19 inline double fluid_density(
const double p,
const double T,
const double x)
26 const double xn = M0 * x / (M0 * x + M1 * (1.0 - x));
29 (M1 * xn + M0 * (1.0 - xn));
38 return 1.0 /
mypow<-i>(x);
41 const double p =
mypow<(i >> 1)>(x);
42 return (i & 1) ? p * p * x : p * p;
53 static double get(
double rho,
double T)
55 const double rho_c = 314;
56 const double CVF = 14.058;
58 const double sigma = 0.36502496e-09;
59 const double k = 1.38062e-23;
60 const double eps = 138.08483e-23;
61 const double c1 = 0.3125;
62 const double c2 = 2.0442e-49;
64 const double T_star = T * k / eps;
67 double Omega = loop1_term<0>(T_star);
68 Omega += loop1_term<1>(T_star);
69 Omega += loop1_term<2>(T_star);
70 Omega += loop1_term<3>(T_star);
71 Omega += loop1_term<4>(T_star);
73 Omega = std::exp(Omega);
76 const double eta_0 = c1 * std::sqrt(c2 * T) / (sigma * sigma * Omega);
78 double sum = loop2_term<2>(rho);
79 sum += loop2_term<3>(rho);
80 sum += loop2_term<4>(rho);
84 CVF * 1e-6 * (
C[0] / (rho -
C[1]) +
C[0] /
C[1] + sum);
93 return A[i] * mypow<i>(log(T_star));
99 return C[i] *
mypow<i - 1>(rho);
102 static const double A[5];
103 static const double C[5];
108 static double get(
double rho,
double T)
110 double my, my_0, my_1;
121 double h[6][7] = {{0.0}};
122 h[0][0] = 0.520094000;
123 h[1][0] = 0.085089500;
124 h[2][0] = -1.083740000;
125 h[3][0] = -0.289555000;
126 h[0][1] = 0.222531000;
127 h[1][1] = 0.999115000;
128 h[2][1] = 1.887970000;
129 h[3][1] = 1.266130000;
130 h[5][1] = 0.120573000;
131 h[0][2] = -0.281378000;
132 h[1][2] = -0.906851000;
133 h[2][2] = -0.772479000;
134 h[3][2] = -0.489837000;
135 h[4][2] = -0.257040000;
136 h[0][3] = 0.161913000;
137 h[1][3] = 0.257399000;
138 h[0][4] = -0.032537200;
139 h[3][4] = 0.069845200;
140 h[4][5] = 0.008721020;
141 h[3][6] = -0.004356730;
142 h[5][6] = -0.000593264;
145 sum1 += H[1] / mypow<1>(T);
146 sum1 += H[2] / mypow<2>(T);
147 sum1 += H[3] / mypow<3>(T);
149 my_0 = 100 * std::sqrt(T) / sum1;
151 double sum2 = inner_loop<0>(rho, T, h);
152 sum2 += inner_loop<1>(rho, T, h);
153 sum2 += inner_loop<2>(rho, T, h);
154 sum2 += inner_loop<3>(rho, T, h);
155 sum2 += inner_loop<4>(rho, T, h);
156 sum2 += inner_loop<5>(rho, T, h);
158 my_1 = std::exp(rho * sum2);
160 my = (my_0 * my_1) / 1e6;
168 const double (&h)[6][7])
170 const double base = rho - 1.0;
172 double sum3 = h[i][0] *
mypow<0>(base);
173 sum3 += h[i][1] * mypow<1>(base);
174 sum3 += h[i][2] * mypow<2>(base);
175 sum3 += h[i][3] * mypow<3>(base);
176 sum3 += h[i][4] * mypow<4>(base);
177 sum3 += h[i][5] * mypow<5>(base);
178 sum3 += h[i][6] * mypow<6>(base);
180 return mypow<i>(1 / T - 1) * sum3;
194 M0 * x / (M0 * x + M1 * (1.0 - x));
197 const double x1 = 1.0 - x0;
200 const double M0_over_M1(M1 / M0);
201 const double V0_over_V1(V0 / V1);
203 const double phi_12 =
205 std::sqrt(V0_over_V1) * std::pow(1.0 / M0_over_M1, 0.25)) /
206 std::sqrt(8.0 * (1.0 + M0_over_M1));
208 const double phi_21 = phi_12 * M0_over_M1 / V0_over_V1;
210 return V0 * x0 / (x0 + x1 * phi_12) + V1 * x1 / (x1 + x0 * phi_21);
215 static double get(
double rho,
double T)
217 const double X1 = 0.95185202;
218 const double X2 = 1.0205422;
220 const double rho_c = 314;
221 const double M = 28.013;
222 const double k = 1.38062e-23;
223 const double eps = 138.08483e-23;
224 const double N_A = 6.02213E26;
225 const double R = 8.31434;
227 const double CCF = 4.173;
229 const double c1 = 0.3125;
230 const double c2 = 2.0442e-49;
231 const double sigma = 0.36502496e-09;
236 const double sum1 = loop1_term<0>(T) + loop1_term<1>(T) +
237 loop1_term<2>(T) + loop1_term<3>(T) +
238 loop1_term<4>(T) + loop1_term<5>(T) +
240 const double temp(std::exp((
f[8] / T)) - 1);
243 (sum1 + ((
f[7] * (
f[8] / T) * (
f[8] / T) * (std::exp((
f[8] / T)))) /
248 cvint = c_v0 * 1000 / N_A;
251 const double log_T_star = std::log(T * k / eps);
254 std::exp(loop2_term<0>(log_T_star) + loop2_term<1>(log_T_star) +
255 loop2_term<2>(log_T_star) + loop2_term<3>(log_T_star) +
256 loop2_term<4>(log_T_star));
260 1e6 * (c1 * std::sqrt(c2 * T) / (sigma * sigma * Omega));
262 const double F = eta_0 * k * N_A / (M * 1000);
264 const double lambda_tr = 2.5 * (1.5 - X1);
265 const double lambda_in = X2 * (cvint / k + X1);
267 const double lambda_0 = F * (lambda_tr + lambda_in);
269 const double sum2 = loop3_term<0>(rho) + loop3_term<1>(rho) +
270 loop3_term<2>(rho) + loop3_term<3>(rho);
271 const double lambda_r = sum2 * CCF;
273 return (lambda_0 + lambda_r) / 1000;
280 return f[i] *
mypow<i - 3>(T);
286 return A[i] * mypow<i>(log_T_star);
292 return C[i] * mypow<i + 1>(rho);
295 const static double A[5];
296 const static double f[9];
297 const static double C[4];
302 static double get(
double rho,
double T)
305 double b[3], B[2], d[4], C[6];
329 const double sum1 = loop_term<0>(T) + loop_term<1>(T) +
330 loop_term<2>(T) + loop_term<3>(T);
332 const double lambda_0 = std::sqrt(T) * sum1;
333 const double lambda_1 =
335 b[2] * std::exp(B[0] * (rho + B[1]) * (rho + B[1]));
337 const double dT = fabs(T - 1) + C[3];
338 const double dT_pow_3_5 = std::pow(dT, 3. / 5.);
339 Q = 2 + (C[4] / dT_pow_3_5);
347 S = C[5] / dT_pow_3_5;
350 const double rho_pow_9_5 = std::pow(rho, 9. / 5.);
351 const double rho_pow_Q = std::pow(rho, Q);
352 const double T_pow_3_2 = T * std::sqrt(T);
353 const double lambda_2 =
354 (d[0] / mypow<10>(T) + d[1]) * rho_pow_9_5 *
355 std::exp(C[0] * (1 - rho * rho_pow_9_5)) +
356 d[2] * S * rho_pow_Q *
357 std::exp((Q / (1. + Q)) * (1 - rho * rho_pow_Q)) +
358 d[3] * std::exp(C[1] * T_pow_3_2 + C[2] / mypow<5>(rho));
360 return lambda_0 + lambda_1 + lambda_2;
364 template <
unsigned i>
367 return a[i] * mypow<i>(T);
370 static const double a[4];
386 const double x0 = std::max(M0 * x / (M0 * x + M1 * (1.0 - x)),
390 const double x1 = 1.0 - x0;
393 const double M1_over_M2 = M1 / M0;
396 const double L1_over_L2 = V1_over_V2 / M1_over_M2;
398 const double M12_pow_mquarter = std::pow(M1_over_M2, -0.25);
399 const double phi_12 = (1.0 + std::sqrt(L1_over_L2) * M12_pow_mquarter) *
400 (1.0 + std::sqrt(V1_over_V2) * M12_pow_mquarter) /
401 std::sqrt(8.0 * (1.0 + M1_over_M2));
402 const double phi_21 = phi_12 * M1_over_M2 / V1_over_V2;
404 return k0 * x0 / (x0 + x1 * phi_12) + k1 * x1 / (x1 + x0 * phi_21);
constexpr double N2
kg mol-1
constexpr double Water
kg mol-1
constexpr double IdealGasConstant
double mypow(const double x)
double fluid_heat_conductivity(const double p, const double T, const double x)
double mypow< 0 >(const double)
double fluid_viscosity(const double p, const double T, const double x)
double fluid_density(const double p, const double T, const double x)
static double loop_term(const double T)
static double get(double rho, double T)
static double loop2_term(const double log_T_star)
static double get(double rho, double T)
static double loop1_term(const double T)
static double loop3_term(const double rho)
static double inner_loop(const double rho, const double T, const double(&h)[6][7])
static double get(double rho, double T)
static double loop2_term(double rho)
static double loop1_term(double T_star)
static double get(double rho, double T)