CoolProp  6.6.1dev
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 #include "CoolPropFluid.h"
4 
5 namespace CoolProp {
6 
8  if (HEOS.is_pure_or_pseudopure) {
9  CoolPropDbl Tstar = HEOS.T() / HEOS.components[0].transport.epsilon_over_k;
10  CoolPropDbl sigma_nm = HEOS.components[0].transport.sigma_eta * 1e9; // 1e9 to convert from m to nm
11  CoolPropDbl molar_mass_kgkmol = HEOS.molar_mass() * 1000; // 1000 to convert from kg/mol to kg/kmol
12 
13  // The nondimensional empirical collision integral from Neufeld
14  // Neufeld, P. D.; Janzen, A. R.; Aziz, R. A. Empirical Equations to Calculate 16 of the Transport Collision Integrals (l,s)*
15  // for the Lennard-Jones (12-6) Potential. J. Chem. Phys. 1972, 57, 1100-1102
16  CoolPropDbl OMEGA22 =
17  1.16145 * pow(Tstar, static_cast<CoolPropDbl>(-0.14874)) + 0.52487 * exp(-0.77320 * Tstar) + 2.16178 * exp(-2.43787 * Tstar);
18 
19  // The dilute gas component -
20  return 26.692e-9 * sqrt(molar_mass_kgkmol * HEOS.T()) / (pow(sigma_nm, 2) * OMEGA22); // Pa-s
21  } else {
22  throw NotImplementedError("TransportRoutines::viscosity_dilute_kinetic_theory is only for pure and pseudo-pure");
23  }
24 }
25 
27  if (HEOS.is_pure_or_pseudopure) {
28  // Retrieve values from the state class
29  CoolProp::ViscosityDiluteGasCollisionIntegralData& data = HEOS.components[0].transport.viscosity_dilute.collision_integral;
30  const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
31  const CoolPropDbl C = data.C, molar_mass = data.molar_mass;
32 
33  CoolPropDbl S;
34  // Unit conversions and variable definitions
35  const CoolPropDbl Tstar = HEOS.T() / HEOS.components[0].transport.epsilon_over_k;
36  const CoolPropDbl sigma_nm = HEOS.components[0].transport.sigma_eta * 1e9; // 1e9 to convert from m to nm
37  const CoolPropDbl molar_mass_kgkmol = molar_mass * 1000; // 1000 to convert from kg/mol to kg/kmol
38 
42  CoolPropDbl summer = 0, lnTstar = log(Tstar);
43  for (std::size_t i = 0; i < a.size(); ++i) {
44  summer += a[i] * pow(lnTstar, t[i]);
45  }
46  S = exp(summer);
47 
48  // The dilute gas component
49  return C * sqrt(molar_mass_kgkmol * HEOS.T()) / (pow(sigma_nm, 2) * S); // Pa-s
50  } else {
51  throw NotImplementedError("TransportRoutines::viscosity_dilute_collision_integral is only for pure and pseudo-pure");
52  }
53 }
54 
56  if (HEOS.is_pure_or_pseudopure) {
57  // Retrieve values from the state class
58  CoolProp::ViscosityDiluteGasPowersOfT& data = HEOS.components[0].transport.viscosity_dilute.powers_of_T;
59  const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
60 
61  CoolPropDbl summer = 0, T = HEOS.T();
62  for (std::size_t i = 0; i < a.size(); ++i) {
63  summer += a[i] * pow(T, t[i]);
64  }
65  return summer;
66  } else {
67  throw NotImplementedError("TransportRoutines::viscosity_dilute_powers_of_T is only for pure and pseudo-pure");
68  }
69 }
71  if (HEOS.is_pure_or_pseudopure) {
72  // Retrieve values from the state class
73  CoolProp::ViscosityDiluteGasPowersOfTr& data = HEOS.components[0].transport.viscosity_dilute.powers_of_Tr;
74  const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
75  CoolPropDbl summer = 0, Tr = HEOS.T() / data.T_reducing;
76  for (std::size_t i = 0; i < a.size(); ++i) {
77  summer += a[i] * pow(Tr, t[i]);
78  }
79  return summer;
80  } else {
81  throw NotImplementedError("TransportRoutines::viscosity_dilute_powers_of_Tr is only for pure and pseudo-pure");
82  }
83 }
84 
86  if (HEOS.is_pure_or_pseudopure) {
87  // Retrieve values from the state class
89  HEOS.components[0].transport.viscosity_dilute.collision_integral_powers_of_Tstar;
90  const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
91 
92  CoolPropDbl summer = 0, Tstar = HEOS.T() / data.T_reducing;
93  for (std::size_t i = 0; i < a.size(); ++i) {
94  summer += a[i] * pow(Tstar, t[i]);
95  }
96  return data.C * sqrt(HEOS.T()) / summer;
97  } else {
98  throw NotImplementedError("TransportRoutines::viscosity_dilute_collision_integral_powers_of_T is only for pure and pseudo-pure");
99  }
100 }
102  if (HEOS.is_pure_or_pseudopure) {
104  HEOS.components[0].transport.viscosity_higher_order.modified_Batschinski_Hildebrand;
105 
106  CoolPropDbl delta = HEOS.rhomolar() / HO.rhomolar_reduce, tau = HO.T_reduce / HEOS.T();
107 
108  // The first term that is formed of powers of tau (Tc/T) and delta (rho/rhoc)
109  CoolPropDbl S = 0;
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]));
112  }
113 
114  // For the terms that multiplies the bracketed term with delta and delta0
115  CoolPropDbl F = 0;
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]);
118  }
119 
120  // for delta_0
121  CoolPropDbl summer_numer = 0;
122  for (unsigned int i = 0; i < HO.g.size(); ++i) {
123  summer_numer += HO.g[i] * pow(tau, HO.h[i]);
124  }
125  CoolPropDbl summer_denom = 0;
126  for (unsigned int i = 0; i < HO.p.size(); ++i) {
127  summer_denom += HO.p[i] * pow(tau, HO.q[i]);
128  }
129  CoolPropDbl delta0 = summer_numer / summer_denom;
130 
131  // The higher-order-term component
132  return S + F * (1 / (delta0 - delta) - 1 / delta0); // Pa-s
133  } else {
134  throw NotImplementedError("TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand is only for pure and pseudo-pure");
135  }
136 }
137 
139  if (HEOS.is_pure_or_pseudopure) {
140  // Retrieve values from the state class
141  CoolProp::ViscosityRainWaterFriendData& data = HEOS.components[0].transport.viscosity_initial.rainwater_friend;
142  const std::vector<CoolPropDbl>&b = data.b, &t = data.t;
143 
144  CoolPropDbl B_eta, B_eta_star;
145  CoolPropDbl Tstar = HEOS.T() / HEOS.components[0].transport.epsilon_over_k; // [no units]
146  CoolPropDbl sigma = HEOS.components[0].transport.sigma_eta; // [m]
147 
148  CoolPropDbl summer = 0;
149  for (unsigned int i = 0; i < b.size(); ++i) {
150  summer += b[i] * pow(Tstar, t[i]);
151  }
152  B_eta_star = summer; // [no units]
153  B_eta = 6.02214129e23 * pow(sigma, 3) * B_eta_star; // [m^3/mol]
154  return B_eta; // [m^3/mol]
155  } else {
156  throw NotImplementedError("TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend is only for pure and pseudo-pure");
157  }
158 }
159 
161  // Inspired by the form from Tariq, JPCRD, 2014
162  if (HEOS.is_pure_or_pseudopure) {
163  // Retrieve values from the state class
164  CoolProp::ViscosityInitialDensityEmpiricalData& data = HEOS.components[0].transport.viscosity_initial.empirical;
165  const std::vector<CoolPropDbl>&n = data.n, &d = data.d, &t = data.t;
166 
167  CoolPropDbl tau = data.T_reducing / HEOS.T(); // [no units]
168  CoolPropDbl delta = HEOS.rhomolar() / data.rhomolar_reducing; // [no units]
169 
170  CoolPropDbl summer = 0;
171  for (unsigned int i = 0; i < n.size(); ++i) {
172  summer += n[i] * pow(delta, d[i]) * pow(tau, t[i]);
173  }
174  return summer; // [Pa-s]
175  } else {
176  throw NotImplementedError("TransportRoutines::viscosity_initial_density_dependence_empirical is only for pure and pseudo-pure");
177  }
178 }
179 
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));
182  double sum;
183  int i, j;
184 
185  // Dilute-gas component
186  *mubar_0 = 100.0 * sqrt(Tbar) / (1.67752 + 2.20462 / Tbar + 0.6366564 / powInt(Tbar, 2) - 0.241605 / powInt(Tbar, 3));
187 
188  //Fill in zeros in H
189  for (i = 0; i <= 5; i++) {
190  for (j = 0; j <= 6; j++) {
191  H[i][j] = 0;
192  }
193  }
194 
195  //Set non-zero parameters of H
196  H[0][0] = 5.20094e-1;
197  H[1][0] = 8.50895e-2;
198  H[2][0] = -1.08374;
199  H[3][0] = -2.89555e-1;
200 
201  H[0][1] = 2.22531e-1;
202  H[1][1] = 9.99115e-1;
203  H[2][1] = 1.88797;
204  H[3][1] = 1.26613;
205  H[5][1] = 1.20573e-1;
206 
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;
212 
213  H[0][3] = 1.61913e-1;
214  H[1][3] = 2.57399e-1;
215 
216  H[0][4] = -3.25372e-2;
217  H[3][4] = 6.98452e-2;
218 
219  H[4][5] = 8.72102e-3;
220 
221  H[3][6] = -4.35673e-3;
222  H[5][6] = -5.93264e-4;
223 
224  // Finite density component
225  sum = 0;
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));
229  }
230  }
231  *mubar_1 = exp(rhobar * sum);
232 }
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));
242  double summer = 0;
243  for (int i = 0; i < 26; ++i) {
244  summer += Bij[i] * pow(1 / Tbar - 1, I[i]) * pow(rhobar - 1, J[i]);
245  }
246  double mu1 = exp(rhobar * summer);
247  double mubar = mu0 * mu1;
248  return 55.2651e-6 * mubar;
249 }
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;
254 
255  pstar = 22.064e6; // [Pa]
256  Tstar = 647.096; // [K]
257  rhostar = 322; // [kg/m^3]
258  Tbar = HEOS.T() / Tstar;
259  rhobar = HEOS.rhomass() / rhostar;
260  R_Water = HEOS.gas_constant() / HEOS.molar_mass(); // [J/kg/K]
261 
262  // Dilute and finite gas portions
263  visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1);
264 
265  // **********************************************************************
266  // ************************ Critical Enhancement ************************
267  // **********************************************************************
268  delta = rhobar;
269  // "Normal" calculation
270  drhodp = 1 / (R_Water * HEOS.T() * (1 + 2 * delta * HEOS.dalphar_dDelta() + delta * delta * HEOS.d2alphar_dDelta2()));
271  drhobar_dpbar = pstar / rhostar * drhodp;
272  // "Reducing" calculation
273  tau = 1 / Tbar_R;
274  drhodp_R = 1
275  / (R_Water * Tbar_R * Tstar
276  * (1 + 2 * rhobar * HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, tau, delta)
277  + delta * delta * HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, tau, delta)));
278  drhobar_dpbar_R = pstar / rhostar * drhodp_R;
279 
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));
285  } else {
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);
288  if (qc * zeta > 1) {
289  L = log((1 + w) / (1 - w));
290  } else {
291  L = 2 * atan(std::abs(w));
292  }
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);
296  }
297  mubar_2 = exp(x_mu * Y);
298 
299  return (mubar_0 * mubar_1 * mubar_2) / 1e6;
300 }
302  CoolPropDbl Tr = HEOS.T() / 591.75, rhor = HEOS.keyed_output(CoolProp::iDmass) / 291.987;
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);
306 }
307 
309  CoolPropDbl Tr = HEOS.T() / 33.145, rhor = HEOS.keyed_output(CoolProp::iDmass) * 0.011;
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));
312 }
314  CoolPropDbl Tr = HEOS.T() / 562.02, rhor = HEOS.rhomass() / 304.792;
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)));
319 }
321 
322  CoolPropDbl Tr = HEOS.T() / 507.82, rhor = HEOS.keyed_output(CoolProp::iDmass) / 233.182;
323 
324  // Output is in Pa-s
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));
330 }
331 
334  CoolPropDbl Tr = HEOS.T() / 540.13, rhor = HEOS.rhomass() / 232;
335 
336  // Output is in Pa-s
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));
341 }
342 
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;
347  // Eq. (9) from Laesecke, JPCRD, 2017
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);
349  // Eq. (8) from Laesecke, JPCRD, 2017
350  double residual = eta_tL * (c1 * Tr * pow(rhor, 3) + (pow(rhor, 2) + pow(rhor, gamma)) / (Tr - c2));
351  return residual;
352 }
353 
355  if (HEOS.is_pure_or_pseudopure) {
356  CoolProp::ViscosityFrictionTheoryData& F = HEOS.components[0].transport.viscosity_higher_order.friction_theory;
357 
358  CoolPropDbl tau = F.T_reduce / HEOS.T(), kii = 0, krrr = 0, kaaa = 0, krr, kdrdr;
359 
360  double psi1 = exp(tau) - F.c1;
361  double psi2 = exp(pow(tau, 2)) - F.c2;
362 
363  double ki = (F.Ai[0] + F.Ai[1] * psi1 + F.Ai[2] * psi2) * tau;
364 
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);
368  if (F.Arr.empty()) {
369  krr = 0;
370  kdrdr = (F.Adrdr[0] + F.Adrdr[1] * psi1 + F.Adrdr[2] * psi2) * pow(tau, F.Nrr);
371  } else {
372  krr = (F.Arr[0] + F.Arr[1] * psi1 + F.Arr[2] * psi2) * pow(tau, F.Nrr);
373  kdrdr = 0;
374  }
375  if (!F.Aii.empty()) {
376  kii = (F.Aii[0] + F.Aii[1] * psi1 + F.Aii[2] * psi2) * pow(tau, F.Nii);
377  }
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);
381  }
382 
383  double p = HEOS.p() / 1e5; // [bar]; 1e5 for conversion from Pa -> bar
384  double pr =
385  HEOS.T() * HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar) / 1e5; // [bar/K]; 1e5 for conversion from Pa -> bar
386  double pa = p - pr; //[bar]
387  double pid = HEOS.rhomolar() * HEOS.gas_constant() * HEOS.T() / 1e5; // [bar]; 1e5 for conversion from Pa -> bar
388  double deltapr = pr - pid;
389 
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;
392 
393  return eta_f; //[Pa-s]
394  } else {
395  throw NotImplementedError("TransportRoutines::viscosity_higher_order_friction_theory is only for pure and pseudo-pure");
396  }
397 }
398 
400  double eta_0, eta_0_slash, eta_E_slash, B, C, D, ln_eta, x;
401  //
402  // Arp, V.D., McCarty, R.D., and Friend, D.G.,
403  // "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
404  // NIST Technical Note 1334 (revised), 1998.
405  //
406  // Using Arp NIST report
407  // Report is not clear on viscosity, referring to REFPROP source code for clarity
408 
409  // Correlation wants density in g/cm^3; kg/m^3 --> g/cm^3, divide by 1000
410  CoolPropDbl rho = HEOS.keyed_output(CoolProp::iDmass) / 1000.0, T = HEOS.T();
411 
412  if (T <= 300) {
413  x = log(T);
414  } else {
415  x = log(300.0);
416  }
417  // Evaluate the terms B,C,D
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;
423 
424  if (T <= 100) {
425  ln_eta = eta_0_slash + eta_E_slash;
426  // 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
427  return exp(ln_eta) / 10.0 / 1e6;
428  } else {
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);
431  // 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
432  return (exp(ln_eta) + eta_0 - exp(eta_0_slash)) / 10.0 / 1e6;
433  }
434 }
435 
437  CoolPropDbl B_eta, C_eta, epsilon_over_k = 577.87, /* [K]*/
438  sigma0 = 0.3408e-9, /* [m] */
439  delta = 0.4575, /* NOT the reduced density, that is rhor here*/
440  N_A = 6.02214129e23, M = 32.04216, /* kg/kmol */
441  T = HEOS.T();
442  CoolPropDbl rhomolar = HEOS.rhomolar();
443 
444  CoolPropDbl B_eta_star, C_eta_star;
445  CoolPropDbl Tstar = T / epsilon_over_k; // [no units]
446  CoolPropDbl rhor = HEOS.rhomass() / 273;
447  CoolPropDbl Tr = T / 512.6;
448 
449  // Rainwater-Friend initial density terms
450  { // Scoped here so that we can re-use the b variable
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};
453  CoolPropDbl summer = 0;
454  for (unsigned int i = 0; i < 9; ++i) {
455  summer += b[i] * pow(Tstar, t[i]);
456  }
457  B_eta_star = summer; // [no units]
458  B_eta = N_A * pow(sigma0, 3) * B_eta_star; // [m^3/mol]
459 
460  CoolPropDbl c[2] = {1.86222085e-3, 9.990338};
461  C_eta_star = c[0] * pow(Tstar, 3) * exp(c[1] * pow(Tstar, static_cast<CoolPropDbl>(-0.5))); // [no units]
462  C_eta = pow(N_A * pow(sigma0, 3), 2) * C_eta_star; // [m^6/mol^2]
463  }
464 
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};
470 
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);
475 
476  CoolPropDbl summerd = 0;
477  for (int i = 0; i < 7; ++i) {
478  summerd += d[i] / pow(Tr, i);
479  }
480  for (int j = 1; j < 10; ++j) {
481  summerd += e[j] * pow(rhor, j);
482  }
483  CoolPropDbl sigmac = 0.7193422e-9; // [m]
484  CoolPropDbl sigma_HS = summerd * sigmac; // [m]
485  CoolPropDbl b = 2 * M_PI * N_A * pow(sigma_HS, 3) / 3; // [m^3/mol]
486  CoolPropDbl zeta = b * rhomolar / 4; // [-]
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); // [-]
489 
490  CoolPropDbl f = 1 / (1 + exp(5 * (rhor - 1)));
491  return eta_0 * (f * eta_g + (1 - f) * eta_E);
492 }
493 
495  double C1 = 1.3163, //
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;
497 
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); // uPa-s
504 
505  double rhobar = HEOS.rhomolar() / 1000; // [mol/L]
506  double eta_L = C2 * (rhoL * rhoL) / (rhoL - rhobar) * sqrt(HEOS.T()) * exp(rhobar / (rhoL - rhobar) * DeltaGstar / (Ru * HEOS.T()));
507 
508  double chi = rhobar - rhocbar;
509  double tau = HEOS.T() - Tc;
510 
511  double DELTAeta_c = 4 * DELTAeta_max / ((exp(chi) + exp(-chi)) * (exp(tau) + exp(-tau)));
512 
513  return (pow((rhoL - rhobar) / rhoL, C1) * eta_DG + pow(rhobar / rhoL, C1) * eta_L + DELTAeta_c) / 1e6;
514 }
515 
517  // From CAO, JPCRD, 2016
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;
523 
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); // [uPa-s]
527 
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; // [uPa-s]
530 
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; // [uPa-s]
534 
535  return (eta0 + eta1 + DELTAeta) / 1e6;
536 }
538  // From CAO, JPCRD, 2016
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;
544 
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); // [uPa-s]
548 
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; // [uPa-s]
551 
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; // [uPa-s]
555 
556  return (eta0 + eta1 + DELTAeta) / 1e6; // [Pa-s]
557 }
559  // From Balogun, JPCRD, 2016
560  double Tr = HEOS.T() / 616.168, rhor = HEOS.rhomolar() / 1000.0 / 2.69392;
561 
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); // [uPa-s]
565 
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; // [uPa-s]
568 
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);
573 
574  return (eta0 + eta1 + DELTAeta) / 1e6; // [Pa-s]
575 }
576 
578  double C[] = {
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;
581 
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);
585  }
586 
587  return 12.0085 * sqrt(Tstar) * OMEGA_2_2 / 1e6; //[Pa-s]
588 }
590  // From Tariq, JPCRD, 2014
591  CoolPropDbl T = HEOS.T();
592  CoolPropDbl S_eta = exp(-1.5093 + 364.87 / T - 39537 / pow(T, 2)); //[nm^2]
593  return 0.19592 * sqrt(T) / S_eta / 1e6; //[Pa-s]
594 }
595 
597  // From Laesecke, JPRCD, 2016
598  double eta0, eta1, DELTAetar, den, Bstar;
599  double T = HEOS.T();
600 
601  double a[] = {
602  1749.354893188350, -369.069300007128, 5423856.34887691, -2.21283852168356, -269503.247933569, 73145.021531826, 5.34368649509278};
603 
604  // Eq. (4) from Laesecke, JPRCD, 2016
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; // [Pa-s]
607  return eta0;
608 }
609 
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};
615 
616  double sum1 = 0, sum2 = 0, tau = 305.33 / HEOS.T(), delta = HEOS.rhomolar() / 6870;
617 
618  for (int i = 1; i <= 9; ++i) {
619  sum1 += g[i] * pow(delta, r[i]) * pow(tau, s[i]);
620  }
621  for (int i = 10; i <= 11; ++i) {
622  sum2 += g[i] * pow(delta, r[i]) * pow(tau, s[i]);
623  }
624  return 15.977 * sum1 / (1 + sum2) / 1e6;
625 }
627  // Retrieve values from the state class
628  CoolProp::ViscosityChungData& data = HEOS.components[0].transport.viscosity_Chung;
629 
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};
634  double A[11];
635 
636  if (HEOS.is_pure_or_pseudopure) {
637  double Vc_cm3mol = 1 / (data.rhomolar_critical / 1e6); // [cm^3/mol]
638  double acentric = data.acentric; // [-]
639  double M_gmol = data.molar_mass * 1000.0; // [g/mol]
640  double Tc = data.T_critical; // [K]
641  double mu_D = data.dipole_moment_D; // [D]
642  double kappa = 0;
643 
644  double mu_r = 131.3 * mu_D / sqrt(Vc_cm3mol * Tc); // [-]
645 
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;
648  }
649  double F_c = 1 - 0.2756 * acentric + 0.059035 * pow(mu_r, 4) + kappa; // [-]
650  double epsilon_over_k = Tc / 1.2593; // [K]
651 
652  double rho_molcm3 = HEOS.rhomolar() / 1e6;
653  double T = HEOS.T();
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; // [P]
658 
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); // [P]
663 
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)); // [P]
666 
667  return (eta_k_P + eta_p_P) / 10.0; // [P] -> [Pa*s]
668  } else {
669  throw NotImplementedError("TransportRoutines::viscosity_Chung is only for pure and pseudo-pure");
670  }
671 }
672 
674  if (HEOS.is_pure_or_pseudopure) {
675  // Retrieve values from the state class
676  CoolProp::ConductivityDiluteRatioPolynomialsData& data = HEOS.components[0].transport.conductivity_dilute.ratio_polynomials;
677 
678  CoolPropDbl summer1 = 0, summer2 = 0, Tr = HEOS.T() / data.T_reducing;
679  for (std::size_t i = 0; i < data.A.size(); ++i) {
680  summer1 += data.A[i] * pow(Tr, data.n[i]);
681  }
682  for (std::size_t i = 0; i < data.B.size(); ++i) {
683  summer2 += data.B[i] * pow(Tr, data.m[i]);
684  }
685 
686  return summer1 / summer2;
687  } else {
688  throw NotImplementedError("TransportRoutines::conductivity_dilute_ratio_polynomials is only for pure and pseudo-pure");
689  }
690 };
691 
693  if (HEOS.is_pure_or_pseudopure) {
694  // Retrieve values from the state class
695  CoolProp::ConductivityResidualPolynomialData& data = HEOS.components[0].transport.conductivity_residual.polynomials;
696 
697  CoolPropDbl summer = 0, tau = data.T_reducing / HEOS.T(), delta = HEOS.keyed_output(CoolProp::iDmass) / data.rhomass_reducing;
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]);
700  }
701  return summer;
702  } else {
703  throw NotImplementedError("TransportRoutines::conductivity_residual_polynomial is only for pure and pseudo-pure");
704  }
705 };
706 
708  if (HEOS.is_pure_or_pseudopure) {
709  // Retrieve values from the state class
711  HEOS.components[0].transport.conductivity_residual.polynomial_and_exponential;
712 
713  CoolPropDbl summer = 0, tau = HEOS.tau(), delta = HEOS.delta();
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]));
716  }
717  return summer;
718  } else {
719  throw NotImplementedError("TransportRoutines::conductivity_residual_polynomial_and_exponential is only for pure and pseudo-pure");
720  }
721 };
722 
724  if (HEOS.is_pure_or_pseudopure) {
725  // Olchowy and Sengers cross-over term
726 
727  // Retrieve values from the state class
728  CoolProp::ConductivityCriticalSimplifiedOlchowySengersData& data = HEOS.components[0].transport.conductivity_critical.Olchowy_Sengers;
729 
730  double k = data.k, R0 = data.R0, nu = data.nu, gamma = data.gamma, GAMMA = data.GAMMA, zeta0 = data.zeta0, qD = data.qD,
731  Tc = HEOS.get_reducing_state().T, // [K]
732  rhoc = HEOS.get_reducing_state().rhomolar, // [mol/m^3]
733  Pcrit = HEOS.get_reducing_state().p, // [Pa]
734  Tref, // [K]
735  cp, cv, delta, num, zeta, mu, pi = M_PI, OMEGA_tilde, OMEGA_tilde0;
736 
737  if (ValidNumber(data.T_ref))
738  Tref = data.T_ref;
739  else
740  Tref = 1.5 * Tc;
741 
742  delta = HEOS.delta();
743 
744  double dp_drho = HEOS.gas_constant() * HEOS.T() * (1 + 2 * delta * HEOS.dalphar_dDelta() + delta * delta * HEOS.d2alphar_dDelta2());
745  double X = Pcrit / pow(rhoc, 2) * HEOS.rhomolar() / dp_drho;
746 
747  double tau_ref = Tc / Tref;
748  double dp_drho_ref = HEOS.gas_constant() * Tref
749  * (1 + 2 * delta * HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, tau_ref, delta)
750  + delta * delta * HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, tau_ref, delta));
751  double Xref = Pcrit / pow(rhoc, 2) * HEOS.rhomolar() / dp_drho_ref * Tref / HEOS.T();
752  num = X - Xref;
753 
754  // No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
755  // See also Lemmon, IJT, 2004, page 27
756  if (num < DBL_EPSILON * 10)
757  return 0.0;
758  else
759  zeta = zeta0 * pow(num / GAMMA, nu / gamma); //[m]
760 
761  cp = HEOS.cpmolar(); //[J/mol/K]
762  cv = HEOS.cvmolar(); //[J/mol/K]
763  mu = HEOS.viscosity(); //[Pa-s]
764 
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))); //[-]
767 
768  double lambda = HEOS.rhomolar() * cp * R0 * k * HEOS.T() / (6 * pi * mu * zeta) * (OMEGA_tilde - OMEGA_tilde0); //[W/m/K]
769  return lambda; //[W/m/K]
770  } else {
771  throw NotImplementedError("TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers is only for pure and pseudo-pure");
772  }
773 };
774 
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));
778 };
779 
781  CoolPropDbl nc = 0.775547504e-3 * 4.81384, Tr = HEOS.T() / 304.1282, alpha, rhor = HEOS.keyed_output(iDmass) / 467.6;
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};
783 
784  // Equation 6 from Scalabrin
785  alpha = 1 - a[10] * acosh(1 + a[11] * pow(pow(1 - Tr, 2), a[12]));
786 
787  // Equation 5 from Scalabrin
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;
792 }
793 
795 
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};
799 
800  //Vesovic Eq. 31 [no units]
801  double summer = 0;
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;
805 
806  //Vesovic Eq. 12 [no units]
807  double r = sqrt(2.0 / 5.0 * cint_k);
808 
809  // According to REFPROP, 1+r^2 = cp-2.5R. This is unclear to me but seems to suggest that cint/k is the difference
810  // between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic
811  // does not yield exactly the correct values
812 
813  Tstar = HEOS.T() / e_k;
814  //Vesovic Eq. 30 [no units]
815  summer = 0;
816  for (int i = 0; i <= 7; i++)
817  summer += b[i] / pow(Tstar, i);
818  double Gstar_lambda = summer;
819 
820  //Vesovic Eq. 29 [W/m/K]
821  double lambda_0 = 0.475598e-3 * sqrt(HEOS.T()) * (1 + r * r) / Gstar_lambda;
822 
823  return lambda_0;
824 }
825 
827 
828  double tau = HEOS.tau();
829  double l[]={0.0151874307, 0.0280674040, 0.0228564190, -0.00741624210};
830  // Huber 2016 Eq. (3)
831  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]
832 
833  return lambda_0/1000;
834 }
835 
837 
838  double e_k = 245.0;
839  double tau = 305.33 / HEOS.T(), Tstar = HEOS.T() / e_k;
840  double fint = 1.7104147 - 0.6936482 / Tstar;
841  double lambda_0 = 0.276505e-3 * (HEOS.calc_viscosity_dilute() * 1e6) * (3.75 - fint * (tau * tau * HEOS.d2alpha0_dTau2() + 1.5)); //[W/m/K]
842 
843  return lambda_0;
844 }
845 
847 
848  if (HEOS.is_pure_or_pseudopure) {
849  CoolProp::ConductivityDiluteEta0AndPolyData& E = HEOS.components[0].transport.conductivity_dilute.eta0_and_poly;
850 
851  double eta0_uPas = HEOS.calc_viscosity_dilute() * 1e6; // [uPa-s]
852  double summer = E.A[0] * eta0_uPas;
853  for (std::size_t i = 1; i < E.A.size(); ++i)
854  summer += E.A[i] * pow(static_cast<CoolPropDbl>(HEOS.tau()), E.t[i]);
855  return summer;
856  } else {
857  throw NotImplementedError("TransportRoutines::conductivity_dilute_eta0_and_poly is only for pure and pseudo-pure");
858  }
859 }
860 
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;
876 }
877 
879 
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}};
885 
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;
888  double xi;
889  int i, j;
890  double R = 461.51805; //[J/kg/K]
891 
892  Tbar = HEOS.T() / Tstar;
893  rhobar = HEOS.keyed_output(CoolProp::iDmass) / rhostar;
894 
895  // Dilute gas contribution
896  lambdabar_0 =
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));
898 
899  sum = 0;
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);
903  }
904  }
905  // Finite density contribution
906  lambdabar_1 = exp(rhobar * sum);
907 
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();
910 
911  double drhodp = 1 / (R * HEOS.T() * (1 + 2 * rhobar * HEOS.dalphar_dDelta() + rhobar * rhobar * HEOS.d2alphar_dDelta2()));
912  double drhobar_dpbar = pstar / rhostar * drhodp;
913  double drhodp_Trbar = 1
914  / (R * Tr_bar * Tstar
915  * (1 + 2 * rhobar * HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, 1 / Tr_bar, delta)
916  + delta * delta * HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, 1 / Tr_bar, delta)));
917  double drhobar_dpbar_Trbar = pstar / rhostar * drhodp_Trbar;
918  double cp = HEOS.cpmass(); // [J/kg/K]
919  double cv = HEOS.cvmass(); // [J/kg/K]
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)
924  xi = 0;
925  else
926  xi = xi_0 * pow(DELTAchibar_T / Lambda_0, nu / gamma);
927  double y = qd_bar * xi;
928 
929  double Z;
930  double kappa = cp / cv;
931  if (y < 1.2e-7)
932  Z = 0;
933  else
934  Z = 2 / (pi * y) * (((1 - 1 / kappa) * atan(y) + y / kappa) - (1 - exp(-1 / (1 / y + y * y / 3 / rhobar / rhobar))));
935 
936  lambdabar_2 = GAMMA * rhobar * cpbar * Tbar / mubar * Z;
937 
938  return (lambdabar_0 * lambdabar_1 + lambdabar_2) * lambdastar;
939 }
940 
942 
943  double B1 = -2.5370, // [mW/m/K]
944  B2 = 0.05366, // [mW/m/K^2]
945  C1 = 0.94215, // [-]
946  C2 = 0.14914, // [mW/m/K^2]
947  DeltaGstar = 2508.58, //[J/mol]
948  rhoL = 68.345, // [mol/dm^3] = [mol/L]
949  rhocbar = 7.5114, // [mol/dm^3]
950  DELTAlambda_max = 25, //[mW/m/K]
951  Ru = 8.31451, // [J/mol/K]
952  Tc = 299.2793, //[K]
953  T = HEOS.T(); //[K]
954 
955  double lambda_DG = B1 + B2 * T;
956 
957  double rhobar = HEOS.rhomolar() / 1000; // [mol/L]
958  double lambda_L = C2 * (rhoL * rhoL) / (rhoL - rhobar) * sqrt(T) * exp(rhobar / (rhoL - rhobar) * DeltaGstar / (Ru * T));
959 
960  double chi = rhobar - rhocbar;
961  double tau = T - Tc;
962 
963  double DELTAlambda_c = 4 * DELTAlambda_max / ((exp(chi) + exp(-chi)) * (exp(tau) + exp(-tau)));
964 
965  return (pow((rhoL - rhobar) / rhoL, C1) * lambda_DG + pow(rhobar / rhoL, C1) * lambda_L + DELTAlambda_c) / 1e3;
966 }
967 
969 
970  /*
971  From "Thermal Conductivity of Ammonia in a Large
972  Temperature and Pressure Range Including the Critical Region"
973  by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
974  Bereicht der Bunsengesellschaft Phys. Chem. 88 (1984) 422-427
975  */
976 
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;
980 
981  rho = HEOS.keyed_output(CoolProp::iDmass);
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; // [Pa-K]
986  X_T = 0.61 * rhoc + 16.5 * log(t);
987  // Along the critical isochore (only a function of temperature) (Eq. 9)
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);
993  } else {
994  DELTA_lambda = DELTA_lambda_id * (X_T * X_T) / (X_T * X_T + powInt(rho - 0.96 * rhoc, 2));
995  }
996 
997  return DELTA_lambda;
998 }
999 
1001  /*
1002  What an incredibly annoying formulation! Implied coefficients?? Not cool.
1003  */
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};
1009  // Equation 17
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);
1013 
1014  // Critical component
1015  lambda_c = 0.0;
1016 
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;
1019 
1020  double DeltaT = std::abs(1 - T / Tc), DeltaRho = std::abs(1 - rho / rhoc_crit);
1021  double eta = HEOS.viscosity(); // [Pa-s]
1022  double K_T = HEOS.isothermal_compressibility(), K_Tprime, K_Tbar;
1024 
1025  double W = pow(DeltaT / 0.2, 2) + pow(DeltaRho / 0.25, 2);
1026 
1027  if (W > 1) {
1028  K_Tbar = K_T;
1029  } else {
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));
1032 
1040  double dhdx =
1041  E1
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));
1045  // Right-hand-side of Equation 9
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;
1049  }
1050 
1051  // 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.
1052  // 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.
1053  lambda_c =
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));
1055  }
1056  return lambda_0 + lambda_e + lambda_c;
1057 }
1059 
1060  double delta = HEOS.rhomolar() / 10139.0, tau = 190.55 / HEOS.T();
1061  double lambda_dilute, lambda_residual, lambda_critical;
1062 
1063  // Viscosity formulation from Friend, JPCRD, 1989
1064  // Dilute
1065  double C[] = {
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);
1071  }
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]);
1080  }
1081  for (int i = 10; i <= 11; ++i) {
1082  summer2 += ge[i] * pow(delta, re[i]) * pow(tau, se[i]);
1083  }
1084  double eta_residual = 12.149 * summer1 / (1 + summer2);
1085  double eta = eta_residual + eta_dilute;
1086 
1087  // 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)); // [mW/m/K]
1090  // Residual
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};
1094  double summer = 0;
1095  for (int i = 1; i <= 6; ++i) {
1096  summer += jl[i] * pow(delta, rl[i]) * pow(tau, sl[i]);
1097  }
1098  double delta_sigma_star = 1.0; // Looks like a typo in Friend - should be 1 instead of 11
1099  if (HEOS.T() < HEOS.T_critical() && HEOS.rhomolar() < HEOS.rhomolar_critical()) {
1100  delta_sigma_star = HEOS.saturation_ancillary(iDmolar, 1, iT, HEOS.T()) / HEOS.keyed_output(CoolProp::irhomolar_critical);
1101  }
1102  lambda_residual = 6.29638 * (summer + jl[7] * POW2(delta) / delta_sigma_star); // [mW/m/K]
1103  // Critical region
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);
1108  double CHI_T_star;
1109  if (std::abs(Tstar) < 0.03) {
1110  if (std::abs(rhostar) < 1e-16) {
1111  // Equation 26
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) {
1115  // Equation 23
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);
1118  double theta = 1;
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);
1121  }
1122  CHI_T_star = Q * pow(std::abs(rhostar), -a) * pow(theta, b) / (theta + OMEGA * (theta + R));
1123  } else {
1124  // Equation 19a
1125  CHI_T_star = 0.28631 * delta * tau / (1 + 2 * delta * HEOS.dalphar_dDelta() + POW2(delta) * HEOS.d2alphar_dDelta2());
1126  }
1127  } else {
1128  // Equation 19a
1129  CHI_T_star = 0.28631 * delta * tau / (1 + 2 * delta * HEOS.dalphar_dDelta() + POW2(delta) * HEOS.d2alphar_dDelta2());
1130  }
1131 
1132  lambda_critical = 91.855 / (eta * POW2(tau)) * POW2(1 + delta * HEOS.dalphar_dDelta() - delta * tau * HEOS.d2alphar_dDelta_dTau())
1133  * pow(CHI_T_star, 0.4681) * F; //[mW/m/K]
1134  return (lambda_dilute + lambda_residual + lambda_critical) * 0.001;
1135 }
1136 
1138  CoolPropDbl& rhomolar0) {
1139  int iter = 0;
1140  double resid = 9e30, resid_old = 9e30;
1141  CoolPropDbl alphar = HEOS.alphar();
1142  CoolPropDbl Z = HEOS.keyed_output(iZ);
1143 
1144  Eigen::Vector2d r;
1145  Eigen::Matrix2d J;
1146  HEOS_Reference.specify_phase(iphase_gas); // Something homogeneous, not checked
1147  // Update the reference fluid with the conformal state
1148  HEOS_Reference.update_DmolarT_direct(rhomolar0, T0);
1149  do {
1150  CoolPropDbl dtau_dT = -HEOS_Reference.T_critical() / (T0 * T0);
1151  CoolPropDbl ddelta_drho = 1 / HEOS_Reference.rhomolar_critical();
1152  // Independent variables are T0 and rhomolar0, residuals are matching alphar and Z
1153  r(0) = HEOS_Reference.alphar() - alphar;
1154  r(1) = HEOS_Reference.keyed_output(iZ) - Z;
1155  J(0, 0) = HEOS_Reference.dalphar_dTau() * dtau_dT;
1156  J(0, 1) = HEOS_Reference.dalphar_dDelta() * ddelta_drho;
1157  // Z = 1+delta*dalphar_ddelta(tau,delta)
1158  // dZ_dT
1159  J(1, 0) = HEOS_Reference.delta() * HEOS_Reference.d2alphar_dDelta_dTau() * dtau_dT;
1160  // dZ_drho
1161  J(1, 1) = (HEOS_Reference.delta() * HEOS_Reference.d2alphar_dDelta2() + HEOS_Reference.dalphar_dDelta()) * ddelta_drho;
1162  // Step in v obtained from Jv = -r
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();
1166  // Calculate the old residual after the last step
1167  resid_old = sqrt(POW2(r(0)) + POW2(r(1)));
1168  for (double frac = 1.0; frac > 0.001; frac /= 2) {
1169  try {
1170  // Calculate new values
1171  double T_new = T0_init + frac * v(0);
1172  double rhomolar_new = rhomolar0_init + frac * v(1);
1173  // Update state with step
1174  HEOS_Reference.update_DmolarT_direct(rhomolar_new, T_new);
1175  resid = sqrt(POW2(HEOS_Reference.alphar() - alphar) + POW2(HEOS_Reference.keyed_output(iZ) - Z));
1176  if (resid > resid_old) {
1177  continue;
1178  }
1179  good_solution = true;
1180  T0 = T_new;
1181  rhomolar0 = rhomolar_new;
1182  break;
1183  } catch (...) {
1184  continue;
1185  }
1186  }
1187  if (!good_solution) {
1188  throw ValueError(format("Not able to get a solution"));
1189  }
1190  iter++;
1191  if (iter > 50) {
1192  throw ValueError(format("conformal_state_solver took too many iterations; residual is %g; prior was %g", resid, resid_old));
1193  }
1194  } while (std::abs(resid) > 1e-9);
1195 }
1196 
1198  // Collect some parameters
1199  CoolPropDbl M = HEOS.molar_mass(), M0 = HEOS_Reference.molar_mass(), Tc = HEOS.T_critical(), Tc0 = HEOS_Reference.T_critical(),
1200  rhocmolar = HEOS.rhomolar_critical(), rhocmolar0 = HEOS_Reference.rhomolar_critical();
1201 
1202  // Get a reference to the ECS data
1203  CoolProp::ViscosityECSVariables& ECS = HEOS.components[0].transport.viscosity_ecs;
1204 
1205  // The correction polynomial psi_eta
1206  double psi = 0;
1207  for (std::size_t i = 0; i < ECS.psi_a.size(); i++) {
1208  psi += ECS.psi_a[i] * pow(HEOS.rhomolar() / ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
1209  }
1210 
1211  // The dilute gas portion for the fluid of interest [Pa-s]
1212  CoolPropDbl eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
1213 
1214  // ************************************
1215  // Start with a guess for theta and phi
1216  // ************************************
1217  CoolPropDbl theta = 1;
1218  CoolPropDbl phi = 1;
1219 
1220  // The equivalent substance reducing ratios
1221  CoolPropDbl f = Tc / Tc0 * theta;
1222  CoolPropDbl h = rhocmolar0 / rhocmolar * phi; // Must be the ratio of MOLAR densities!!
1223 
1224  // To be solved for
1225  CoolPropDbl T0 = HEOS.T() / f;
1226  CoolPropDbl rhomolar0 = HEOS.rhomolar() * h;
1227 
1228  // **************************
1229  // Solver for conformal state
1230  // **************************
1231 
1232  //
1233  HEOS_Reference.specify_phase(iphase_gas); // something homogeneous
1234 
1235  conformal_state_solver(HEOS, HEOS_Reference, T0, rhomolar0);
1236 
1237  // Update the reference fluid with the updated conformal state
1238  HEOS_Reference.update_DmolarT_direct(rhomolar0 * psi, T0);
1239 
1240  // Recalculate ESRR
1241  f = HEOS.T() / T0;
1242  h = rhomolar0 / HEOS.rhomolar(); // Must be the ratio of MOLAR densities!!
1243 
1244  // **********************
1245  // Remaining calculations
1246  // **********************
1247 
1248  // The reference fluid's contribution to the viscosity [Pa-s]
1249  CoolPropDbl eta_resid = HEOS_Reference.calc_viscosity_background();
1250 
1251  // The F factor
1252  CoolPropDbl F_eta = sqrt(f) * pow(h, -static_cast<CoolPropDbl>(2.0L / 3.0L)) * sqrt(M / M0);
1253 
1254  // The total viscosity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
1255  CoolPropDbl eta = eta_dilute + eta_resid * F_eta;
1256 
1257  return eta;
1258 }
1259 
1261 
1262  // Get a reference to the data
1263  const CoolProp::ViscosityRhoSrVariables& data = HEOS.components[0].transport.viscosity_rhosr;
1264 
1265  // The dilute gas portion for the fluid of interest [Pa-s]
1266  CoolPropDbl eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
1267 
1268  // Calculate x
1269  double x = HEOS.rhomolar() * HEOS.gas_constant() * (HEOS.tau() * HEOS.dalphar_dTau() - HEOS.alphar()) / data.rhosr_critical;
1270 
1271  // Crossover variable
1272  double psi_liq = 1 / (1 + exp(-100.0 * (x - 2)));
1273 
1274  // Evaluated using Horner's method
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])));
1278 
1279  // Evaluate the reference fluid
1280  double etastar_ref = exp(psi_liq * f_liq + (1 - psi_liq) * f_vap);
1281 
1282  // Get the non-dimensionalized viscosity
1283  double etastar_fluid = 1 + data.C * (etastar_ref - 1);
1284 
1285  return etastar_fluid * eta_dilute;
1286 }
1287 
1289  // Collect some parameters
1290  CoolPropDbl M = HEOS.molar_mass(), M_kmol = M * 1000, M0 = HEOS_Reference.molar_mass(), Tc = HEOS.T_critical(), Tc0 = HEOS_Reference.T_critical(),
1291  rhocmolar = HEOS.rhomolar_critical(), rhocmolar0 = HEOS_Reference.rhomolar_critical(), R_u = HEOS.gas_constant(),
1292  R = HEOS.gas_constant() / HEOS.molar_mass(), //[J/kg/K]
1293  R_kJkgK = R_u / M_kmol;
1294 
1295  // Get a reference to the ECS data
1296  CoolProp::ConductivityECSVariables& ECS = HEOS.components[0].transport.conductivity_ecs;
1297 
1298  // The correction polynomial psi_eta in rho/rho_red
1299  double psi = 0;
1300  for (std::size_t i = 0; i < ECS.psi_a.size(); ++i) {
1301  psi += ECS.psi_a[i] * pow(HEOS.rhomolar() / ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
1302  }
1303 
1304  // The correction polynomial f_int in T/T_red
1305  double fint = 0;
1306  for (std::size_t i = 0; i < ECS.f_int_a.size(); ++i) {
1307  fint += ECS.f_int_a[i] * pow(HEOS.T() / ECS.f_int_T_reducing, ECS.f_int_t[i]);
1308  }
1309 
1310  // The dilute gas density for the fluid of interest [uPa-s]
1311  CoolPropDbl eta_dilute = viscosity_dilute_kinetic_theory(HEOS) * 1e6;
1312 
1313  // The mass specific ideal gas constant-pressure specific heat [J/kg/K]
1314  CoolPropDbl cp0 = HEOS.calc_cpmolar_idealgas() / HEOS.molar_mass();
1315 
1316  // The internal contribution to the thermal conductivity [W/m/K]
1317  CoolPropDbl lambda_int = fint * eta_dilute * (cp0 - 2.5 * R) / 1e3;
1318 
1319  // The dilute gas contribution to the thermal conductivity [W/m/K]
1320  CoolPropDbl lambda_dilute = 15.0e-3 / 4.0 * R_kJkgK * eta_dilute;
1321 
1322  // ************************************
1323  // Start with a guess for theta and phi
1324  // ************************************
1325 
1326  CoolPropDbl theta = 1;
1327  CoolPropDbl phi = 1;
1328 
1329  // The equivalent substance reducing ratios
1330  CoolPropDbl f = Tc / Tc0 * theta;
1331  CoolPropDbl h = rhocmolar0 / rhocmolar * phi; // Must be the ratio of MOLAR densities!!
1332 
1333  // Initial values for the conformal state
1334  CoolPropDbl T0 = HEOS.T() / f;
1335  CoolPropDbl rhomolar0 = HEOS.rhomolar() * h;
1336 
1337  // **************************
1338  // Solver for conformal state
1339  // **************************
1340 
1341  try {
1342  conformal_state_solver(HEOS, HEOS_Reference, T0, rhomolar0);
1343  } catch (std::exception& e) {
1344  throw ValueError(format("Conformal state solver failed; error was %s", e.what()));
1345  }
1346 
1347  // Update the reference fluid with the conformal state
1348  HEOS_Reference.update(DmolarT_INPUTS, rhomolar0 * psi, T0);
1349 
1350  // Recalculate ESRR
1351  f = HEOS.T() / T0;
1352  h = rhomolar0 / HEOS.rhomolar(); // Must be the ratio of MOLAR densities!!
1353 
1354  // The reference fluid's contribution to the conductivity [W/m/K]
1355  CoolPropDbl lambda_resid = HEOS_Reference.calc_conductivity_background();
1356 
1357  // The F factor
1358  CoolPropDbl F_lambda = sqrt(f) * pow(h, static_cast<CoolPropDbl>(-2.0 / 3.0)) * sqrt(M0 / M);
1359 
1360  // The critical contribution from the fluid of interest [W/m/K]
1362 
1363  // The total thermal conductivity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
1364  CoolPropDbl lambda = lambda_int + lambda_dilute + lambda_resid * F_lambda + lambda_critical;
1365 
1366  return lambda;
1367 }
1368 
1369 }; /* namespace CoolProp */