Cantera  3.2.0a2
Loading...
Searching...
No Matches
HighPressureGasTransport.cpp
Go to the documentation of this file.
1/**
2 * @file HighPressureGasTransport.cpp
3 * Implementation file for class HighPressureGasTransport
4 **/
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
15#include <boost/algorithm/string.hpp>
16
17namespace Cantera
18{
19
20/**
21 * @brief Returns interpolated value of (DP)_R obtained from the data in Table 2 of
22 * the Takahashi 1975 paper, given a value of the reduced pressure (Pr) and reduced
23 * temperature (Tr).
24 *
25 * @param Pr Reduced pressure
26 * @param Tr Reduced temperature
27 */
28double takahashiCorrectionFactor(double Pr, double Tr)
29{
30 // In the low pressure limit, no correction is needed. Interpolate
31 // the value towards 1 as pressure drops below the 0.1 threshold.
32 if (Pr < 0.1) {
33 return 1.0;
34 }
35
36 // Data from Table 2 of Takahashi 1975 paper:
37 const static double Pr_lookup[17] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0,
38 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0};
39 const static double DP_Rt_lookup[17] = {1.01, 1.01, 1.01, 1.01, 1.01, 1.01,
40 1.01, 1.02, 1.02, 1.02, 1.02, 1.03, 1.03, 1.04, 1.05, 1.06, 1.07};
41 const static double A_ij_lookup[17] = {0.038042, 0.067433, 0.098317,
42 0.137610, 0.175081, 0.216376, 0.314051, 0.385736, 0.514553, 0.599184,
43 0.557725, 0.593007, 0.696001, 0.790770, 0.502100, 0.837452, 0.890390};
44 const static double B_ij_lookup[17] = {1.52267, 2.16794, 2.42910, 2.77605,
45 2.98256, 3.11384, 3.50264, 3.07773, 3.54744, 3.61216, 3.41882, 3.18415,
46 3.37660, 3.27984, 3.39031, 3.23513, 3.13001};
47 const static double C_ij_lookup[17] = {0., 0., 0., 0., 0., 0., 0., 0.141211,
48 0.278407, 0.372683, 0.504894, 0.678469, 0.665702, 0., 0.602907, 0., 0.};
49 const static double E_ij_lookup[17] = {1., 1., 1., 1., 1., 1., 1., 13.45454,
50 14., 10.00900, 8.57519, 10.37483, 11.21674, 1., 6.19043, 1., 1.};
51
52 // Interpolate to obtain the value of (DP)_R at
53 // the provided value of the reduced pressure (Pr).
54 int Pr_lower = 0; // Index of the lower bounding value of Pr
55 int Pr_upper = 0; // Index of the upper bounding value of Pr
56 double frac = 0.0;
57
58 bool found = false;
59 for (int j = 1; j < 17; j++) {
60 if (Pr_lookup[j] > Pr) {
61 frac = (Pr - Pr_lookup[j-1]) / (Pr_lookup[j] - Pr_lookup[j-1]);
62 found = true;
63 Pr_lower = j-1;
64 Pr_upper = j;
65 break;
66 }
67 }
68 // If this loop completes without finding a bounding value of Pr, use
69 // the final table value.
70 if (!found) {
71 Pr_lower = 16;
72 Pr_upper = 16;
73 frac = 1.0;
74 }
75
76 // Compute the value of (DP)_R at the given Pr value by interpolating the
77 // bounding values of (DP)_R.
78 double DP_Rt = DP_Rt_lookup[Pr_lower];
79 double A = A_ij_lookup[Pr_lower];
80 double B = B_ij_lookup[Pr_lower];
81 double C = C_ij_lookup[Pr_lower];
82 double E = E_ij_lookup[Pr_lower];
83
84 double DP_R_lower = DP_Rt*(1.0 - A*pow(Tr,-B))*(1.0 - C*pow(Tr,-E));
85
86 DP_Rt = DP_Rt_lookup[Pr_upper];
87 A = A_ij_lookup[Pr_upper];
88 B = B_ij_lookup[Pr_upper];
89 C = C_ij_lookup[Pr_upper];
90 E = E_ij_lookup[Pr_upper];
91
92 double DP_R_upper = DP_Rt*(1.0 - A*pow(Tr,-B))*(1.0 - C*pow(Tr,-E));
93
94 // Linear interpolation of the two bounding values of (DP)_R.
95 return DP_R_lower*(1.0 - frac) + DP_R_upper*frac;
96}
97
99{
100 // Call the base class's method to fill the properties
102
103 // Contents of 'critical-properties.yaml', loaded later if needed
104 AnyMap critPropsDb;
105 std::unordered_map<string, AnyMap*> dbSpecies;
106
107 // If a species has a zero acentric factor, check the critical-properties.yaml
108 // database to see if it has a value specified for the species.
109 for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
110 if (m_w_ac[k] == 0.0) {
111 // Load 'crit-properties.yaml' file if not already loaded
112 if (critPropsDb.empty()) {
113 critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
114 dbSpecies = critPropsDb["species"].asMap("name");
115 }
116
117 // All names in critical-properties.yaml are upper case
118 auto ucName = boost::algorithm::to_upper_copy(m_thermo->species(k)->name);
119 if (dbSpecies.count(ucName)) {
120 auto& spec = *dbSpecies.at(ucName);
121 auto& critProps = spec["critical-parameters"].as<AnyMap>();
122 if (critProps.hasKey("acentric-factor")) {
123 m_w_ac[k] = critProps.convert("acentric-factor", "1");
124 }
125 }
126 }
127 }
128}
129
131{
132 size_t nSpecies = m_thermo->nSpecies();
133 m_Tcrit.resize(nSpecies);
134 m_Pcrit.resize(nSpecies);
135 m_Vcrit.resize(nSpecies);
136 m_Zcrit.resize(nSpecies);
137
138 vector<double> molefracs(nSpecies);
139 m_thermo->getMoleFractions(&molefracs[0]);
140
141 vector<double> mf_temp(nSpecies, 0.0);
142
143 for (size_t i = 0; i < nSpecies; ++i) {
144 mf_temp[i] = 1.0;
145 m_thermo->setMoleFractions(&mf_temp[0]);
146
147 if (m_thermo->critTemperature() > 1e4) {
148 throw CanteraError(
149 "HighPressureGasTransport::initializeCriticalProperties",
150 "Species '{}' must have critical properties defined or non-zero "
151 "cubic parameters. Check the species definition or the thermo "
152 "data file.",
153 m_thermo->species(i)->name);
154 }
159
160 mf_temp[i] = 0.0; // Reset for the next iteration
161 }
162
163 // Restore actual mole fractions
164 m_thermo->setMoleFractions(&molefracs[0]);
165}
166
167// Pure species critical properties - Tc, Pc, Vc, Zc:
169{
170 return m_Tcrit[i];
171}
172
174{
175 return m_Pcrit[i];
176}
177
179{
180 return m_Vcrit[i];
181}
182
184{
185 return m_Zcrit[i];
186}
187
189 for (size_t i = 0; i < m_nsp; i++) {
190 for (size_t j = 0; j < m_nsp; j++) {
191 // Add an offset to avoid a condition where x_i and x_j both equal
192 // zero (this would lead to Pr_ij = Inf).
193 double x_i = std::max(Tiny, m_molefracs[i]);
194 double x_j = std::max(Tiny, m_molefracs[j]);
195
196 // Weight mole fractions of i and j so that X_i + X_j = 1.0.
197 double sum_x_ij = x_i + x_j;
198 x_i = x_i/(sum_x_ij);
199 x_j = x_j/(sum_x_ij);
200
201 // Calculate Tr and Pr based on mole-fraction-weighted critical constants.
202 double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
203 double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
204
205 // Calculate the parameters for Takahashi correlation
206 double P_corr_ij;
207 P_corr_ij = takahashiCorrectionFactor(Pr_ij, Tr_ij);
208
209 // If the reduced temperature is too low, the correction factor
210 // P_corr_ij will be < 0.
211 if (P_corr_ij<0) {
212 P_corr_ij = Tiny;
213 }
214 m_P_corr_ij(i, j) = P_corr_ij;
215 }
216 }
217}
218
219void HighPressureGasTransportBase::getBinaryDiffCoeffs(const size_t ld, double* const d)
220{
221 update_C();
222 update_T();
224 // If necessary, evaluate the binary diffusion coefficients from the polynomial
225 // fits
226 if (!m_bindiff_ok) {
227 updateDiff_T();
228 }
229 if (ld < m_nsp) {
230 throw CanteraError("HighPressureGasTransport::getBinaryDiffCoeffs",
231 "ld is too small");
232 }
233
234 double rp = 1.0/m_thermo->pressure();
235 for (size_t i = 0; i < m_nsp; i++) {
236 for (size_t j = 0; j < m_nsp; j++) {
237 // Multiply the standard low-pressure binary diffusion coefficient
238 // (m_bdiff) by the Takahashi correction factor P_corr_ij.
239 d[ld*j + i] = m_P_corr_ij(i,j)*(rp * m_bdiff(i,j));
240 }
241 }
242}
243
245{
246 update_T();
247 update_C();
249
250 // update the binary diffusion coefficients if necessary
251 if (!m_bindiff_ok) {
252 updateDiff_T();
253 }
254
255 double mmw = m_thermo->meanMolecularWeight();
256 double p = m_thermo->pressure();
257 if (m_nsp == 1) {
258 d[0] = m_P_corr_ij(0,0)*m_bdiff(0,0) / p;
259 } else {
260 for (size_t i = 0; i < m_nsp; i++) {
261 double sum2 = 0.0;
262 for (size_t j = 0; j < m_nsp; j++) {
263 if (j != i) {
264 sum2 += m_molefracs[j] / (m_P_corr_ij(i,j)*m_bdiff(j,i));
265 }
266 }
267 if (sum2 <= 0.0) {
268 d[i] = m_P_corr_ij(i,i)*m_bdiff(i,i) / p;
269 } else {
270 d[i] = (mmw - m_molefracs[i] * m_mw[i])/(p * mmw * sum2);
271 }
272 }
273 }
274}
275
277{
278 update_T();
279 update_C();
281
282 // update the binary diffusion coefficients if necessary
283 if (!m_bindiff_ok) {
284 updateDiff_T();
285 }
286
287 double p = m_thermo->pressure();
288 if (m_nsp == 1) {
289 d[0] = m_P_corr_ij(0,0)*m_bdiff(0,0) / p;
290 } else {
291 for (size_t i = 0; i < m_nsp; i++) {
292 double sum2 = 0.0;
293 for (size_t j = 0; j < m_nsp; j++) {
294 if (j != i) {
295 sum2 += m_molefracs[j] / (m_P_corr_ij(i,j)*m_bdiff(j,i));
296 }
297 }
298 if (sum2 <= 0.0) {
299 d[i] = m_P_corr_ij(i,i)*m_bdiff(i,i) / p;
300 } else {
301 d[i] = (1 - m_molefracs[i]) / (p * sum2);
302 }
303 }
304 }
305}
306
308{
309 update_T();
310 update_C();
312
313 // update the binary diffusion coefficients if necessary
314 if (!m_bindiff_ok) {
315 updateDiff_T();
316 }
317
318 double mmw = m_thermo->meanMolecularWeight();
319 double p = m_thermo->pressure();
320
321 if (m_nsp == 1) {
322 d[0] = m_P_corr_ij(0,0)*m_bdiff(0,0) / p;
323 } else {
324 for (size_t i=0; i<m_nsp; i++) {
325 double sum1 = 0.0;
326 double sum2 = 0.0;
327 for (size_t j=0; j<m_nsp; j++) {
328 if (j==i) {
329 continue;
330 }
331 sum1 += m_molefracs[j] / (m_P_corr_ij(i,j)*m_bdiff(i,j));
332 sum2 += m_molefracs[j] * m_mw[j] / (m_P_corr_ij(i,j)*m_bdiff(i,j));
333 }
334 sum1 *= p;
335 sum2 *= p * m_molefracs[i] / (mmw - m_mw[i]*m_molefracs[i]);
336 d[i] = 1.0 / (sum1 + sum2);
337 }
338 }
339}
340
341// HighPressureGasTransport Implementation
342// ---------------------------------------
344{
348}
349
351{
352 update_T();
353 vector<double> molefracs(m_nsp);
354 m_thermo->getMoleFractions(&molefracs[0]);
355 vector<double> cp_0_R(m_nsp);
356 m_thermo->getCp_R_ref(&cp_0_R[0]); // Cp/R
357
358 // A model constant from the Euken correlation for polyatomic gases, described
359 // below Equation 1 in ely-hanley1981 .
360 const double f_int = 1.32;
361
362 // Pure-species model parameters
363 // Internal contribution to thermal conductivity (lamba'')
364 vector<double> Lambda_1_i(m_nsp);
365 vector<double> f_i(m_nsp);
366 vector<double> h_i(m_nsp);
367 vector<double> V_k(m_nsp);
368
369 m_thermo -> getPartialMolarVolumes(&V_k[0]);
370 for (size_t i = 0; i < m_nsp; i++) {
371 // Calculate variables for density-independent component, Equation 1,
372 // the equation requires the pure-species viscosity estimate from
373 // Ely and Hanley.
374 double mu_i = elyHanleyDilutePureSpeciesViscosity(V_k[i], Tcrit_i(i),
375 Vcrit_i(i), Zcrit_i(i),
376 m_w_ac[i], m_mw[i]);
377
378 // This is the internal contribution to the thermal conductivity of
379 // pure-species component, i, from Equation 1 in ely-hanley1983
380 Lambda_1_i[i] = (mu_i / m_mw[i])*f_int*GasConstant*(cp_0_R[i] - 2.5);
381
382 // Calculate variables for density-dependent component (lambda')
383 double Tr = m_thermo->temperature() / Tcrit_i(i);
384 double Vr = V_k[i] / Vcrit_i(i);
385 double theta_i = thetaShapeFactor(Tr, Vr, m_w_ac[i]);
386 double phi_i = phiShapeFactor(Tr, Vr, Zcrit_i(i), m_w_ac[i]);
387
388 f_i[i] = (Tcrit_i(i) / m_ref_Tc)*theta_i; // Equation 12 ely-hanley1983
389 h_i[i] = (Vcrit_i(i) / m_ref_Vc)*phi_i; // Equation 13 ely-hanley1983
390 }
391
392 double h_m = 0; // Corresponding states parameter, h_x,0 from ely-hanley1983
393 double f_m = 0; // Corresponding states parameter, f_x,0 from ely-hanley1983
394 double mw_m = 0; // Mixture molecular weight
395 double Lambda_1_m = 0.0; // Internal component of mixture thermal conductivity
396 for (size_t i = 0; i < m_nsp; i++) {
397 for (size_t j = 0; j < m_nsp; j++) {
398 // Compute the internal contribution to the thermal conductivity of the
399 // mixture
400
401 // Equation 3 ely-hanley1983
402 double Lambda_1_ij = 2*Lambda_1_i[i]*Lambda_1_i[j] /
403 (Lambda_1_i[i] + Lambda_1_i[j] + Tiny);
404 // Equation 2, ely-hanley1983
405 Lambda_1_m += molefracs[i]*molefracs[j]*Lambda_1_ij;
406
407 // Variables for density-dependent translational/collisional component of
408 // the mixture.
409
410 // Equation 10, ely-hanley1983
411 double f_ij = sqrt(f_i[i]*f_i[j]);
412
413 // Equation 11, ely-hanley1983
414 double h_ij = 0.125*pow(pow(h_i[i],1.0/3.0) + pow(h_i[j],1.0/3.0), 3.0);
415
416 // Equation 15, ely-hanley1983
417 double mw_ij_inv = 0.5*(m_mw[i] + m_mw[j])/(m_mw[i]*m_mw[j]);
418
419 // Equation 8, ely-hanley1983
420 f_m += molefracs[i]*molefracs[j]*f_ij*h_ij;
421
422 // Equation 9, ely-hanley1983
423 h_m += molefracs[i]*molefracs[j]*h_ij;
424
425 // Equation, 14 ely-hanley1983
426 mw_m += molefracs[i]*molefracs[j]*sqrt(mw_ij_inv*f_ij)*pow(h_ij,-4.0/3.0);
427 }
428 }
429
430 // The following two equations are the final steps for Equations 8 and 14 in
431 // ely-hanley1983 . The calculations in the loop above computed the
432 // right-hand-side of the equations, but the left-hand-side of the equations
433 // contain other variables that must be moved to the right-hand-side in order to
434 // get the values of the variables of interest.
435 f_m = f_m/h_m;
436 mw_m = pow(mw_m,-2.0)*f_m*pow(h_m,-8.0/3.0);
437
438 // The two relations below are from Equation 7, ely-hanley1983 . This must
439 // be in units of g/cm^3 for use with the empirical correlation.
440 const double kg_m3_to_g_cm3 = 1e-3; // Conversion factor from kg/m^3 to g/cm^3
441 double rho_0 = m_thermo->density()*h_m*kg_m3_to_g_cm3;
442 double T_0 = m_temp/f_m;
443
444 // Equation 18, ely-hanley1983
445 double Lambda_2_ref = elyHanleyReferenceThermalConductivity(rho_0, T_0);
446
447 // Equation 6, ely-hanley1983
448 double F_m = sqrt(f_m*m_ref_MW/mw_m)*pow(h_m,-2.0/3.0);
449
450 // Equation 5, ely-hanley1983
451 double Lambda_2_m = F_m*Lambda_2_ref;
452
453 return Lambda_1_m + Lambda_2_m;
454}
455
457 double Tc, double Vc, double Zc, double acentric_factor, double mw)
458{
459 double Tr = m_thermo->temperature() / Tc;
460 double Vr = V / Vc;
461 double theta_i = thetaShapeFactor(Tr, Vr, acentric_factor);
462 double phi_i = phiShapeFactor(Tr, Vr, Zc, acentric_factor);
463
464 double f_fac = (Tc / m_ref_Tc)*theta_i; // Equation 7 ely-hanley1981
465 double h_fac = (Vc / m_ref_Vc)*phi_i; // Equation 8 ely-hanley1981
466 double T_0 = m_temp/f_fac; // Equation 3, ely-hanley1981
467
468 // Dilute reference fluid viscosity correlation, from Table III in
469 // ely-hanley1981
470 double mu_0 = elyHanleyDiluteReferenceViscosity(T_0);
471
472 // Equation 2, ely-hanley1981
473 double F = sqrt(f_fac*(mw/m_ref_MW))*pow(h_fac,-2.0/3.0);
474
475 return mu_0*F;
476}
477
479 double acentric_factor)
480{
481 double T_p = std::min(std::max(Tr,0.5), 2.0);
482 double V_p = std::min(std::max(Vr,0.5), 2.0);
483
484 return 1 + (acentric_factor - m_ref_acentric_factor)*(0.090569 - 0.862762*log(T_p)
485 + (0.316636 - 0.465684/T_p)*(V_p - 0.5));
486
487}
488
489double HighPressureGasTransport::phiShapeFactor(double Tr, double Vr, double Zc,
490 double acentric_factor)
491{
492 double T_p = std::min(std::max(Tr,0.5), 2.0);
493 double V_p = std::min(std::max(Vr,0.5), 2.0);
494
495 return (1 + (acentric_factor - m_ref_acentric_factor)*(0.394901*(V_p - 1.023545)
496 - 0.932813*(V_p - 0.754639)*log(T_p)))*(m_ref_Zc/Zc);
497
498}
499
501{
502 // Conversion factor from the correlation in micrograms/cm/s to Pa*s.
503 double const correlation_viscosity_conversion = 1e-7;
504
505 if (T0 > 10000) {
506 T0 = 10000; // Limit the temperature to 10000 K
507 }
508
509 // Coefficients for the correlation from Table III of ely-hanley1981
510 const vector<double> c = {2.907741307e6, -3.312874033e6, 1.608101838e6,
511 -4.331904871e5, 7.062481330e4, -7.116620750e3,
512 4.325174400e2, -1.445911210e1, 2.037119479e-1};
513
514 double mu_0 = 0.0;
515 for (size_t i = 0; i < 9; i++) {
516 mu_0 += c[i]*pow(T0,(i+1.0-4.0)/3.0);
517 }
518 return correlation_viscosity_conversion*mu_0;
519}
520
522 double T0)
523{
524 // Computing the individual terms of Equation 18, ely-hanley1983 . This is an
525 // evaluation of the expressions shown in Table I of ely-hanley1983 . The
526 // correlation returns values of thermal conductivity in mW/m/K,
527 // so a conversion is needed.
528 const double correlation_factor = 1e-3; // mW/m/K to W/m/K
529
530 // This is the reference gas, dilute gas viscosity (eta_0 in Table III of
531 // ely-hanley1981)
532 double mu_0 = elyHanleyDiluteReferenceViscosity(T0);
533
534 // First term in Equation 18. This expression has the correct units because
535 // it does not use any empirical correlation, so it is excluded at the end from
536 // the unit conversion.
537 double Lambda_ref_star = (15*GasConstant / (4*m_ref_MW))*mu_0;
538
539 // Second term in Equation 18
540 const vector<double> b = {-2.52762920e-1, 3.34328590e-1, 1.12, 1.680e2};
541 double Lambda_ref_1 = (b[0] + b[1]*pow(b[2] - log(T0/b[3]), 2))*rho0;
542
543 // Third term in Equation 18
544 const vector<double> a = {-7.197708227, 8.5678222640e1, 1.2471834689e1,
545 -9.8462522975e2, 3.5946850007e-1, 6.9798412538e1,
546 -8.7288332851e2};
547 double delta_lambda_ref = exp(a[0] + a[1]/T0)
548 * (exp((a[2] + a[3]*pow(T0,-1.5))*pow(rho0,0.1)
549 + (rho0/m_ref_rhoc - 1)*sqrt(rho0)*(a[4] + a[5]/T0
550 + a[6]*pow(T0,-2))) - 1.0);
551
552 return Lambda_ref_star + (Lambda_ref_1 + delta_lambda_ref)*correlation_factor;
553}
554
556{
558
559 // This is η*ξ
560 double nondimensional_viscosity = highPressureNondimensionalViscosity(
563
564 // Using equation 9-4.14, with units of 1/(Pa*s)
565 double numerator = GasConstant*m_Tc_mix*pow(Avogadro,2.0);
566 double denominator = pow(m_MW_mix,3.0)*pow(m_Pc_mix,4.0);
567 double xi = pow(numerator / denominator, 1.0/6.0);
568
569 // Return the viscosity in kg/m/s
570 return nondimensional_viscosity / xi;
571}
572
574{
575 double Tc_mix = 0.0;
576 double Pc_mix_n = 0.0; // Numerator in equation 9-5.19
577 double Pc_mix_d = 0.0; // Denominator in equation 9-5.19
578
579 // Equation 9.5.20, Cantera already mole-weights the molecular weights
580 double MW_mix = m_thermo->meanMolecularWeight();
581
582 // Mole-fraction-weighted mixture average of the low-pressure polarity correction
583 // factor
584 double FP_mix_o = 0;
585
586 // Mole-fraction-weighted mixture average of the low-pressure quantum correction
587 // factor
588 double FQ_mix_o = 0;
589
590 double MW_H = m_mw[0]; // Molecular weight of the heaviest species
591 double MW_L = m_mw[0]; // Molecular weight of the lightest species
592
593 double tKelvin = m_thermo->temperature();
594 double P_vap_mix = m_thermo->satPressure(tKelvin);
595 size_t nsp = m_thermo->nSpecies();
596 vector<double> molefracs(nsp);
597 m_thermo->getMoleFractions(&molefracs[0]);
598
599 double x_H = molefracs[0]; // Holds the mole fraction of the heaviest species
600 for (size_t i = 0; i < m_nsp; i++) {
601 // Calculate pure-species critical constants and add their contribution
602 // to the mole-fraction-weighted mixture averages:
603 double Tc = Tcrit_i(i);
604 double Tr = tKelvin/Tc;
605 double Zc = Zcrit_i(i);
606
607 Tc_mix += Tc*molefracs[i]; // Equation 9-5.18
608
609 Pc_mix_n += molefracs[i]*Zc; // Numerator of 9-5.19
610 // (Units of Vcrit_i are m^3/kmol, which are fine because they cancel out
611 // with the Cantera gas constant's units used later)
612 Pc_mix_d += molefracs[i]*Vcrit_i(i); // Denominator of 9-5.19
613
614 // Calculate ratio of heaviest to lightest species, used in
615 // Equation 9-5.23.
616 if (m_mw[i] > MW_H) {
617 MW_H = m_mw[i];
618 x_H = molefracs[i];
619 } else if (m_mw[i] < MW_L) {
620 MW_L = m_mw[i];
621 }
622
623 // Calculate pure-species reduced dipole moment for pure-species
624 // polar correction term.
625 // Equation 9-4.17 requires the pressure
626 // to be in units of bar, so we convert from Pa to bar.
627 // The dipole moment is stored in SI units, and it needs to be in
628 // units of Debye for the Lucas method.
629 double pascals_to_bar = 1e-5;
630 double SI_to_Debye = lightSpeed / 1e-21; // Conversion factor from C*m to Debye
631 double dipole_ii = m_dipole(i,i)*SI_to_Debye;
632 double mu_ri = 52.46*dipole_ii*dipole_ii*(Pcrit_i(i)*pascals_to_bar)/(Tc*Tc);
633
634 // mole-fraction weighting of pure-species polar correction term
635 FP_mix_o += molefracs[i] * polarityCorrectionFactor(mu_ri, Tr, Zc);
636
637 // Calculate contribution to quantum correction term.
638 // Note: This assumes the species of interest (He, H2, and D2) have
639 // been named in this specific way.
640 vector<string> spnames = m_thermo->speciesNames();
641 if (spnames[i] == "He") {
642 FQ_mix_o += molefracs[i]*quantumCorrectionFactor(1.38, Tr, m_mw[i]);
643 } else if (spnames[i] == "H2") {
644 FQ_mix_o += molefracs[i]*(quantumCorrectionFactor(0.76, Tr, m_mw[i]));
645 } else if (spnames[i] == "D2") {
646 FQ_mix_o += molefracs[i]*(quantumCorrectionFactor(0.52, Tr, m_mw[i]));
647 } else {
648 FQ_mix_o += molefracs[i];
649 }
650 }
651
652 double Tr_mix = tKelvin/Tc_mix;
653 double Pc_mix = GasConstant*Tc_mix*Pc_mix_n/Pc_mix_d;
654 double Pr_mix = m_thermo->pressure()/Pc_mix;
655
656 // Compute the mixture value of the low-pressure quantum correction factor
657 // Equation 9-5.23.
658 double ratio = MW_H/MW_L;
659 double A = 1.0;
660 if (ratio > 9 && x_H > 0.05 && x_H < 0.7) {
661 A = 1 - 0.01*pow(ratio,0.87);
662 }
663 FQ_mix_o *= A;
664
665 m_FQ_mix_o = FQ_mix_o;
666 m_FP_mix_o = FP_mix_o;
667 m_Tc_mix = Tc_mix;
668 m_Tr_mix = Tr_mix;
669 m_Pc_mix = Pc_mix;
670 m_Pr_mix = Pr_mix;
671 m_MW_mix = MW_mix;
672 m_P_vap_mix = P_vap_mix;
673}
674
676 double Tr, double FP, double FQ)
677{
678 double first_term = 0.807*pow(Tr,0.618) - 0.357*exp(-0.449*Tr);
679 double second_term = 0.340*exp(-4.058*Tr) + 0.018;
680 return (first_term + second_term)*FP*FQ;
681}
682
684 double Tr, double Pr, double FP_low, double FQ_low, double P_vap, double P_crit)
685{
686 // This is η_0*ξ
687 double Z_1 = lowPressureNondimensionalViscosity(Tr, FP_low, FQ_low);
688
689 double Z_2;
690 if (Tr <= 1.0) {
691 if (Pr < P_vap/P_crit) {
692 double alpha = 3.262 + 14.98*pow(Pr, 5.508);
693 double beta = 1.390 + 5.746*Pr;
694 Z_2 = 0.600 + 0.760*pow(Pr,alpha) + (0.6990*pow(Pr,beta) - 0.60) * (1-Tr);
695 } else {
696 throw CanteraError(
697 "HighPressureGasTransport::highPressureNondimensionalViscosity",
698 "State is outside the limits of the Lucas model, Pr ({}) >= "
699 "P_vap / P_crit ({}) when Tr ({}) <= 1.0", Pr, P_vap / P_crit, Tr);
700 }
701 } else if (Tr > 1.0 && Tr < 40.0) {
702 if (Pr > 0.0 && Pr <= 100.0) {
703 // The following expressions are given in page 9.36 of Poling and
704 // correspond to parameters in equation 9-6.8.
705 double a_1 = 1.245e-3;
706 double a_2 = 5.1726;
707 double gamma = -0.3286;
708 double a = a_1*exp(a_2*pow(Tr,gamma))/Tr;
709
710 double b_1 = 1.6553;
711 double b_2 = 1.2723;
712 double b = a*(b_1*Tr - b_2);
713
714 double c_1 = 0.4489;
715 double c_2 = 3.0578;
716 double delta = -37.7332;
717 double c = c_1*exp(c_2*pow(Tr, delta))/Tr;
718
719 double d_1 = 1.7368;
720 double d_2 = 2.2310;
721 double epsilon = -7.6351;
722 double d = d_1*exp(d_2*pow(Tr, epsilon))/Tr;
723
724 double e = 1.3088;
725
726 double f_1 = 0.9425;
727 double f_2 = -0.1853;
728 double zeta = 0.4489;
729 double f = f_1*exp(f_2*pow(Tr, zeta));
730
731 Z_2 = Z_1*(1 + (a*pow(Pr,e)) / (b*pow(Pr,f) + pow(1+c*pow(Pr,d),-1)));
732 } else {
733 throw CanteraError(
734 "HighPressureGasTransport::highPressureNondimensionalViscosity",
735 "The value of Pr ({}) is outside the limits of the Lucas model, "
736 "valid values of Pr are: 0.0 < Pr <= 100", Pr);
737 }
738 } else {
739 throw CanteraError(
740 "HighPressureGasTransport::highPressureNondimensionalViscosity",
741 "The value of Tr is outside the limits of the Lucas model, "
742 "valid values of Tr are: 1.0 < Tr < 40", Tr);
743 }
744
745 double Y = Z_2 / Z_1;
746 double FP = (1 + (FP_low - 1)*pow(Y,-3.0)) / FP_low;
747 double FQ = (1 + (FQ_low - 1)*(1.0/Y - 0.007*pow(log(Y),4.0))) / FQ_low;
748
749 // Return the non-dimensional viscosity η*ξ
750 return Z_2 * FP * FQ;
751}
752
754 double MW)
755{
756 return 1.22*pow(Q,0.15)*(1 + 0.00385*pow(pow(Tr - 12.0, 2.0), 1.0/MW)
757 *sign(Tr - 12.0));
758}
759
761 double Z_crit)
762{
763 if (mu_r < 0.022) {
764 return 1;
765 } else if (mu_r < 0.075) {
766 return 1 + 30.55*pow(std::max(0.292 - Z_crit,0.0), 1.72);
767 } else {
768 return 1 + 30.55*pow(std::max(0.292 - Z_crit, 0.0), 1.72)
769 *fabs(0.96 + 0.1*(Tr - 0.7));
770 }
771}
772
773
774// ChungHighPressureGasTransport Implementation
775// --------------------------------------------
777{
782}
783
785{
786 // First fill the species-specific values that will then be used in the
787 // combining rules for the Chung method.
788 m_sigma_i.resize(m_nsp);
791 m_MW_i.resize(m_nsp);
792 m_kappa_i.resize(m_nsp);
793 for (size_t i = 0; i < m_nsp; i++) {
794 // From equation 9-5.32.
795 double m3_per_kmol_to_cm3_per_mol = 1e3; // Convert from m^3/kmol to cm^3/mol
796 double Vc = Vcrit_i(i) * m3_per_kmol_to_cm3_per_mol;
797 m_sigma_i[i] = 0.809*pow(Vc, 1.0/3.0);
798
799 // From equation 9-5.34.
800 m_epsilon_over_k_i[i] = Tcrit_i(i)/1.2593;
801
802 // NOTE: The association parameter is assumed to be zero for all species, but
803 // is left here for completeness or future revision.
804 m_kappa_i[i] = 0.0;
805
806 // These values are available from the base class
808 m_MW_i[i] = m_mw[i];
809 }
810}
811
813{
815
816 // Compute T_star using equation 9-5.26, using the mixture parameters
817 double tKelvin = m_thermo->temperature();
818 double T_star = tKelvin / m_epsilon_over_k_mix;
819
820 // The density is required for high-pressure gases.
821 // The Chung method requires density to be units of mol/cm^3
822 // Use the mixture molecular weight (units of kg/kmol).
823 // 1 kmol/m^3 = 1e-3 mol/cm^3
824 double kg_per_m3_to_mol_per_cm3 = (1.0 / m_MW_mix)*1e-3;
825 double density = m_thermo->density()*kg_per_m3_to_mol_per_cm3;
826
827 // The value of Cv is already a mole-weighted average of the pure species values
828 double Cv_mix = m_thermo->cv_mole(); // Units are J/kmol/K
829
830 // This result is in units of W/m/K
831 double thermal_conductivity = highPressureThermalConductivity(
832 tKelvin, T_star, m_MW_mix, density, Cv_mix, m_Vc_mix,
835
836 // Return the thermal conductivity in W/m/K
837 return thermal_conductivity;
838}
839
841 double T, double T_star, double MW, double rho, double Cv, double Vc,
842 double Tc, double sigma, double acentric_factor, double mu_r,
843 double kappa)
844{
845 // Calculate the low-pressure viscosity using the Chung method (units of Pa*s)
846 double viscosity = lowPressureViscosity(T, T_star, MW, acentric_factor, mu_r,
847 sigma, kappa);
848
849
850 double M_prime = MW / 1000.0; // Convert kg/kmol to kg/mol
851
852 // Definition of tabulated coefficients for the Chung method, as shown in
853 // Table 10-3 on page 10.23.
854 static const vector<double> a = {2.44166, -5.0924e-1, 6.6107, 1.4543e1, 7.9274e-1,
855 -5.8634, 9.1089e1};
856 static const vector<double> b = {7.4824e-1, -1.5094, 5.6207, -8.9139, 8.2019e-1,
857 1.2801e1, 1.2811e2};
858 static const vector<double> c = {-9.1858e-1, -4.9991e1, 6.4760e1, -5.6379,
859 -6.9369e-1, 9.5893, -5.4217e1};
860 static const vector<double> d = {1.2172e2, 6.9983e1, 2.7039e1, 7.4344e1, 6.3173,
861 6.5529e1, 5.2381e2};
862
863 // This is slightly pedantic, but this is done to have the naming convention in the
864 // equations used match the variable names in the code. This is equation 10-5.9.
865 double B_1 = a[0] + b[0]*acentric_factor + c[0]*pow(mu_r, 4.0) + d[0]*kappa;
866 double B_2 = a[1] + b[1]*acentric_factor + c[1]*pow(mu_r, 4.0) + d[1]*kappa;
867 double B_3 = a[2] + b[2]*acentric_factor + c[2]*pow(mu_r, 4.0) + d[2]*kappa;
868 double B_4 = a[3] + b[3]*acentric_factor + c[3]*pow(mu_r, 4.0) + d[3]*kappa;
869 double B_5 = a[4] + b[4]*acentric_factor + c[4]*pow(mu_r, 4.0) + d[4]*kappa;
870 double B_6 = a[5] + b[5]*acentric_factor + c[5]*pow(mu_r, 4.0) + d[5]*kappa;
871 double B_7 = a[6] + b[6]*acentric_factor + c[6]*pow(mu_r, 4.0) + d[6]*kappa;
872
873 double y = rho*Vc/6.0; // Equation 10-5.6 (with rho = 1/V)
874
875 double G_1 = (1.0 - 0.5*y)/(pow(1.0-y, 3.0)); // Equation 10-5.7
876
877 // Equation 10-5.8
878 double G_2 = (B_1*((1.0-exp(-B_4*y)) / y) + B_2*G_1*exp(B_5*y) + B_3*G_1)
879 / (B_1*B_4 + B_2 + B_3);
880
881 double q = 3.586e-3*sqrt(Tc/M_prime) / pow(Vc, 2.0/3.0); // Below equation 10-5.5
882
883 double Tr = T/Tc; // Reduced temperature
884 // The following 4 equations are defined below Equation 10-3.14
885 double alpha = (Cv / GasConstant) - 3.0/2.0; // GasConstant(R) has units (J/kmol/K)
886 double beta = 0.7862 - 0.7109*acentric_factor
887 + 1.3168*acentric_factor*acentric_factor;
888 double Z = 2.0 + 10.5*Tr*Tr;
889 double psi = 1.0 + alpha*((0.215 + 0.28288*alpha - 1.061*beta + 0.26665*Z)
890 / (0.6366 + beta*Z + 1.061*alpha*beta));
891
892 // Equation 10-5.5
893 double lambda = (31.2*viscosity*psi/M_prime)*(1.0/G_2 + B_6*y)
894 + q*B_7*y*y*sqrt(Tr)*G_2;
895
896 // Units are W/m/K
897 return lambda;
898}
899
901{
903
904 // Compute T_star using equation 9-5.26, using the mixture parameters
905 double tKelvin = m_thermo->temperature();
906 double T_star = tKelvin / m_epsilon_over_k_mix;
907
908 // The density is required for high-pressure gases.
909 // The Chung method requires density to be units of mol/cm^3
910 // Use the mixture molecular weight (units of kg/kmol) here.
911 // 1 kmol/m^3 = 1e-3 mol/cm^3
912 double kg_per_m3_to_mol_per_cm3 = (1.0 / m_MW_mix)*1e-3;
913 double molar_density = m_thermo->density()*kg_per_m3_to_mol_per_cm3;
914
915 // This result is in units of micropoise
916 double viscosity = highPressureViscosity(T_star, m_MW_mix, molar_density,
920
921 double micropoise_to_pascals_second = 1e-7;
922 return viscosity*micropoise_to_pascals_second;
923}
924
926{
927 // Here we use the combining rules defined on page 9.25.
928 // We have ASSUMED that the binary interaction parameters are unity for all species
929 // as was done in the Chung method.
930 //
931 // The sigma & kappa relations can be fully computed in the loop. The other ones
932 // require a final division by the final mixture values, and so the quantities in
933 // the loop are only the numerators of the equations that are referenced in the
934 // comments. After the loop, the final mixture values are computed and stored.
935
936 // Zero out the mixture values of the parameters
937 m_sigma_mix = 0.0;
939 m_MW_mix = 0.0;
941 m_mu_mix = 0.0;
942 m_kappa_mix = 0.0;
943 m_Tc_mix = 0.0;
944 m_Vc_mix = 0.0;
945 m_mu_r_mix = 0.0;
946
947 vector<double> molefracs(m_nsp);
948 m_thermo->getMoleFractions(&molefracs[0]);
949 for (size_t i = 0; i < m_nsp; i++) {
950 for (size_t j = 0; j <m_nsp; j++){
951 double sigma_ij = sqrt(m_sigma_i[i]*m_sigma_i[j]); // Equation 9-5.33
952 // Equation 9-5.25
953 m_sigma_mix += molefracs[i]*molefracs[j]*pow(sigma_ij,3.0);
954
955 // Equation 9-5.35
956 double epsilon_over_k_ij = sqrt(m_epsilon_over_k_i[i]*m_epsilon_over_k_i[j]);
957
958 // Equation 9-5.27 (numerator only)
959 m_epsilon_over_k_mix += molefracs[i]*molefracs[j]*epsilon_over_k_ij
960 *pow(sigma_ij,3.0);
961
962 // This equation is raised to the power 2, so what we do here is first
963 // store only the double summation into this variable, and then later
964 // square the result.
965 // Equation 9-5.40
966 double MW_ij = (2*m_MW_i[i]*m_MW_i[j]) / (m_MW_i[i] + m_MW_i[j]);
967
968 // Equation 9-5.28 (double summation only)
969 m_MW_mix += molefracs[i]*molefracs[j]*epsilon_over_k_ij*pow(sigma_ij,2.0)
970 *sqrt(MW_ij);
971
972 // Equation 9-5.36
973 double acentric_factor_ij = 0.5*(m_acentric_factor_i[i]
975
976 // Equation 9-5.29
977 m_acentric_factor_mix += molefracs[i]*molefracs[j]*acentric_factor_ij
978 *pow(sigma_ij,3.0);
979
980 // The base class' dipole moment values are in the SI units, so we need to
981 // convert to Debye units.
982 double SI_to_Debye = lightSpeed / 1e-21; // Conversion from C*m to Debye
983 double dipole_ii = m_dipole(i,i)*SI_to_Debye;
984 double dipole_jj = m_dipole(j,j)*SI_to_Debye;
985
986 // Equation 9-5.30
987 m_mu_mix += molefracs[i]*molefracs[j]*pow(dipole_ii*dipole_jj,2.0)
988 /pow(sigma_ij,3.0);
989
990 // Using equation 9-5.31
991 double kappa_ij = sqrt(m_kappa_i[i]*m_kappa_i[j]); // Equation 9-5.39
992 m_kappa_mix += molefracs[i]*molefracs[j]*kappa_ij; // Equation 9-5.31
993 }
994 }
995
996 // Finalize the expressions for the mixture values of the parameters
997
998 // Equation 9-5.25 computed the cube of sigma_mix
999 m_sigma_mix = pow(m_sigma_mix, 1.0/3.0);
1001
1002 // The MW_mix was only the numerator inside the brackets of equation 9-5.28
1004
1006
1007 // Equation 9-5.30 computed the 4th power of mu_mix
1008 m_mu_mix = pow(pow(m_sigma_mix, 3.0)*m_mu_mix, 1.0/4.0);
1009
1010 // Tc_mix is computed using equation 9-5.44.
1012
1013 // Vc_mix is computed using equation 9-5.43.
1014 m_Vc_mix = pow(m_sigma_mix/0.809, 3.0);
1015
1016 // mu_r_mix is computed using equation 9-5.42.
1017 m_mu_r_mix = 131.3*m_mu_mix/sqrt(m_Vc_mix*m_Tc_mix);
1018}
1019
1020/**
1021 * Returns the value of the Neufeld collision integral for a given dimensionless
1022 * temperature. Implementation of equation 9-4.3.
1023 * Applicable over the range of 0.3 <= T_star <= 100.
1024 *
1025 * @param T_star Dimensionless temperature (Defined in Equation 9-4.1)
1026 */
1027double neufeldCollisionIntegral(double T_star)
1028{
1029 double A = 1.16145;
1030 double B = 0.14874;
1031 double C = 0.52487;
1032 double D = 0.77320;
1033 double E = 2.16178;
1034 double F = 2.43787;
1035
1036 return A / pow(T_star, B) + C / exp(D*T_star) + E / exp(F*T_star);
1037}
1038
1040 double MW, double acentric_factor, double mu_r, double sigma, double kappa)
1041{
1042 double omega = neufeldCollisionIntegral(T_star); // Equation 9-4.3.
1043
1044 // Molecular shapes and polarities factor, Equation 9-4.11
1045 double Fc = 1 - 0.2756*acentric_factor + 0.059035*pow(mu_r, 4.0) + kappa;
1046
1047 // Equation 9-3.9, multiplied by the Chung factor, Fc
1048 // (another way of writing 9-4.10 that avoids explicit use of the critical volume
1049 // in this method)
1050 double viscosity = Fc*(26.69*sqrt(MW*T)/(sigma*sigma*omega));
1051
1052 double micropoise_to_pascals_second = 1e-7;
1053 return micropoise_to_pascals_second*viscosity;
1054}
1055
1057 double rho, double Vc, double Tc, double acentric_factor, double mu_r,
1058 double kappa)
1059{
1060 // Definition of tabulated coefficients for the Chung method, as shown in Table 9-6
1061 // on page 9.40 of Poling.
1062 static const vector<double> a = {6.324, 1.210e-3, 5.283, 6.623, 19.745, -1.900,
1063 24.275, 0.7972, -0.2382, 0.06863};
1064 static const vector<double> b = {50.412, -1.154e-3, 254.209, 38.096, 7.630,
1065 -12.537, 3.450, 1.117, 0.06770, 0.3479};
1066 static const vector<double> c = {-51.680, -6.257e-3, -168.48, -8.464, -14.354,
1067 4.958, -11.291, 0.01235, -0.8163, 0.5926};
1068 static const vector<double> d ={1189.0, 0.03728, 3898.0, 31.42, 31.53, -18.15,
1069 69.35, -4.117, 4.025, -0.727};
1070
1071 // This is slightly pedantic, but this is done to have the naming convention in the
1072 // equations used match the variable names in the code.
1073 double E_1 = a[0] + b[0]*acentric_factor + c[0]*pow(mu_r, 4.0) + d[0]*kappa;
1074 double E_2 = a[1] + b[1]*acentric_factor + c[1]*pow(mu_r, 4.0) + d[1]*kappa;
1075 double E_3 = a[2] + b[2]*acentric_factor + c[2]*pow(mu_r, 4.0) + d[2]*kappa;
1076 double E_4 = a[3] + b[3]*acentric_factor + c[3]*pow(mu_r, 4.0) + d[3]*kappa;
1077 double E_5 = a[4] + b[4]*acentric_factor + c[4]*pow(mu_r, 4.0) + d[4]*kappa;
1078 double E_6 = a[5] + b[5]*acentric_factor + c[5]*pow(mu_r, 4.0) + d[5]*kappa;
1079 double E_7 = a[6] + b[6]*acentric_factor + c[6]*pow(mu_r, 4.0) + d[6]*kappa;
1080 double E_8 = a[7] + b[7]*acentric_factor + c[7]*pow(mu_r, 4.0) + d[7]*kappa;
1081 double E_9 = a[8] + b[8]*acentric_factor + c[8]*pow(mu_r, 4.0) + d[8]*kappa;
1082 double E_10 = a[9] + b[9]*acentric_factor + c[9]*pow(mu_r, 4.0) + d[9]*kappa;
1083
1084 double y = rho*Vc/6.0; // Equation 9-6.20
1085
1086 double G_1 = (1.0 - 0.5*y)/(pow(1.0-y, 3.0)); // Equation 9-6.21
1087
1088 // Equation 9-6.22
1089 double G_2 = (E_1*((1.0-exp(-E_4*y)) / y) + E_2*G_1*exp(E_5*y) + E_3*G_1)
1090 / (E_1*E_4 + E_2 + E_3);
1091
1092 // Equation 9-6.23
1093 double eta_2 = E_7*y*y*G_2*exp(E_8 + E_9/T_star + E_10/(T_star*T_star));
1094
1095 double omega = neufeldCollisionIntegral(T_star); // Equation 9-4.3
1096
1097 // Molecular shapes and polarities factor, Equation 9-4.11
1098 double Fc = 1 -0.2756*acentric_factor + 0.059035*pow(mu_r, 4.0) + kappa;
1099
1100 // Renamed eta_star and eta_star_star from the Poling description to eta_1 and
1101 // eta_2 for naming simplicity.
1102 // Equation 9-6.19
1103 double eta_1 = (sqrt(T_star)/omega) * (Fc*(1.0/G_2 + E_6*y)) + eta_2;
1104
1105 double eta = eta_1 * 36.344 * sqrt(MW*Tc) / pow(Vc, 2.0/3.0); // Equation 9-6.18
1106
1107 return eta;
1108}
1109
1110}
Interface for class HighPressureGasTransport.
ThermoPhase object for the ideal gas equation of state - workhorse for Cantera (see Thermodynamic Pro...
Declaration for class Cantera::Species.
Header file defining class TransportFactory (see TransportFactory)
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition AnyMap.cpp:1468
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition AnyMap.cpp:1841
Base class for exceptions thrown by Cantera classes.
double m_kappa_mix
Mixture association factor [unitless].
vector< double > m_MW_i
Molecular weight [kg/kmol].
void initializePureFluidProperties()
Computes and stores pure-fluid specific properties each species.
double thermalConductivity() override
Calculates the high-pressure mixture thermal conductivity using the Chung method.
double m_MW_mix
Effective mixture molecular weight [kg/kmol].
double m_mu_r_mix
Mixture reduced dipole moment [unitless].
double m_Tc_mix
Mixture critical temperature [K].
double m_epsilon_over_k_mix
Mixture characteristic temperature [K].
vector< double > m_epsilon_over_k_i
Characteristic temperature [K].
vector< double > m_acentric_factor_i
Acentric factor [unitless].
double highPressureThermalConductivity(double T, double T_star, double MW, double rho, double Cv, double Vc, double Tc, double sigma, double acentric_factor, double mu_r, double kappa)
Computes the high-pressure thermal conductivity [W/m/K] using the Chung method.
void init(ThermoPhase *thermo, int mode=0) override
Initialize a transport manager.
double lowPressureViscosity(double T, double T_star, double MW, double acentric_factor, double mu_r, double sigma, double kappa)
Returns the low-pressure mixture viscosity [Pa·s] using the Chung method.
double m_Vc_mix
Mixture critical volume [m³/kmol].
vector< double > m_sigma_i
Effective molecular diameter [Angstroms].
double m_mu_mix
Mixture dipole moment [Debye].
double viscosity() override
Returns the high-pressure mixture viscosity [Pa·s] using the Chung method.
double m_sigma_mix
Effective mixture molecular diameter [Angstroms].
double m_acentric_factor_mix
Mixture acentric factor [unitless].
vector< double > m_kappa_i
Association factor [unitless].
double highPressureViscosity(double T_star, double MW, double rho, double Vc, double Tc, double acentric_factor, double mu_r, double kappa)
Returns the high-pressure mixture viscosity [micropoise] using the Chung method.
void computeMixtureParameters()
Computes the composition-dependent values of the parameters needed for the Chung viscosity model.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
virtual void getTransportData()
Read the transport database.
vector< double > m_mw
Local copy of the species molecular weights.
vector< double > m_molefracs
Vector of species mole fractions.
double m_temp
Current value of the temperature [K] at which the properties in this object are calculated.
virtual void updateDiff_T()
Update the binary diffusion coefficients.
bool m_bindiff_ok
Update boolean for the binary diffusivities at unit pressure.
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
DenseMatrix m_dipole
The effective dipole moment [Coulomb·m] for (i,j) collisions.
vector< double > m_w_ac
Pitzer acentric factor [dimensionless].
void getBinaryDiffCoeffs(const size_t ld, double *const d) override
Computes the matrix of binary diffusion coefficients [m²/s] using the Takahashi correction factor.
vector< double > m_Pcrit
Critical pressure [Pa] of each species.
DenseMatrix m_P_corr_ij
Matrix of Takahashi binary diffusion coefficient corrections.
virtual void updateCorrectionFactors()
Updates the matrix of species-pair Takahashi correction factors for use in computing the binary diffu...
double Pcrit_i(size_t i)
Returns the stored value of the critical pressure for species 'i'.
void getMixDiffCoeffs(double *const d) override
Returns the mixture-averaged diffusion coefficients [m²/s].
double Vcrit_i(size_t i)
Returns the stored value of the critical volume for species 'i'.
void getMixDiffCoeffsMole(double *const d) override
Returns the mixture-averaged diffusion coefficients [m²/s].
vector< double > m_Zcrit
Critical compressibility [unitless] of each species.
double Tcrit_i(size_t i)
Returns the stored value of the critical temperature for species 'i'.
double Zcrit_i(size_t i)
Returns the stored value of the critical compressibility for species 'i'.
vector< double > m_Vcrit
Critical volume [m³/kmol] of each species.
void getTransportData() override
Obtain required parameters from the 'critical-parameters' species input section, and checks the criti...
vector< double > m_Tcrit
Critical temperature [K] of each species.
void getMixDiffCoeffsMass(double *const d) override
Returns the mixture-averaged diffusion coefficients [m²/s].
void initializeCriticalProperties()
Computes and stores the estimate of the critical properties for each species.
double polarityCorrectionFactor(double mu_r, double Tr, double Z_c)
Returns the polarity correction term for a species based on reduced temperature, reduced dipole momen...
const double m_ref_Vc
Critical volume [m^3/kmol].
double elyHanleyDilutePureSpeciesViscosity(double V, double Tc, double Vc, double Zc, double acentric_factor, double mw)
Get the viscosity [Pa·s] of a pure species using the method of Ely and Hanley.
double elyHanleyDiluteReferenceViscosity(double T0)
Returns the viscosity [Pa·s] for the reference fluid (methane) for low pressures.
double thermalConductivity() override
Returns the mixture high-pressure thermal conductivity [W/m/K] using a method of corresponding states...
double highPressureNondimensionalViscosity(double Tr, double Pr, double FP_low, double FQ_low, double P_vap, double P_crit)
Returns the non-dimensional high-pressure mixture viscosity in using the Lucas method.
double quantumCorrectionFactor(double Q, double Tr, double MW)
Calculates quantum correction term of the Lucas method for a species based on the reduced temperature...
const double m_ref_MW
Molecular weight [kg/kmol].
double lowPressureNondimensionalViscosity(double Tr, double FP, double FQ)
Returns the non-dimensional low-pressure mixture viscosity in using the Lucas method.
double elyHanleyReferenceThermalConductivity(double rho0, double T0)
Returns the thermal conductivity [W/m/K] of the reference fluid of methane.
double m_FP_mix_o
Polarity correction factor.
void init(ThermoPhase *thermo, int mode=0) override
Initialize a transport manager.
const double m_ref_Zc
Critical compressibility [unitless].
double viscosity() override
Returns the mixture high-pressure viscosity [Pa·s] using the Lucas method.
void computeMixtureParameters()
Computes the composition-dependent values of parameters that are needed for the Lucas viscosity model...
double m_FQ_mix_o
Quantum correction factor.
const double m_ref_Tc
Critical temperature [K].
double phiShapeFactor(double Tr, double Vr, double Zc, double acentric_factor)
Returns the phi shape factor of Leach and Leland for a pure species.
const double m_ref_acentric_factor
Acentric factor [unitless].
double thetaShapeFactor(double Tr, double Vr, double acentric_factor)
Returns the theta shape factor of Leach and Leland for a pure species.
const double m_ref_rhoc
Critical density [g/cm^3].
void update_T() override
Update the internal parameters whenever the temperature has changed.
void init(ThermoPhase *thermo, int mode=0) override
Initialize a transport manager.
void update_C() override
Update the internal parameters whenever the concentrations have changed.
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
Definition Phase.cpp:289
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:232
double temperature() const
Temperature (K).
Definition Phase.h:563
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:656
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:434
virtual double density() const
Density (kg/m^3).
Definition Phase.h:588
const vector< string > & speciesNames() const
Return a const reference to the vector of species names.
Definition Phase.cpp:148
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:874
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:581
Base class for a phase with thermodynamic properties.
virtual double critTemperature() const
Critical temperature (K).
virtual void getCp_R_ref(double *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
virtual double critPressure() const
Critical pressure (Pa).
virtual double critVolume() const
Critical volume (m3/kmol).
virtual double cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
virtual double satPressure(double t)
Return the saturation pressure given the temperature.
virtual double critCompressibility() const
Critical compressibility (unitless).
ThermoPhase * m_thermo
pointer to the object representing the phase
Definition Transport.h:426
size_t m_nsp
Number of species in the phase.
Definition Transport.h:429
ThermoPhase & thermo()
Phase object.
Definition Transport.h:103
int sign(T x)
Sign of a number.
Definition global.h:339
const double Avogadro
Avogadro's Number [number/kmol].
Definition ct_defs.h:81
const double lightSpeed
Speed of Light in a vacuum [m/s].
Definition ct_defs.h:93
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const double Tiny
Small number to compare differences of mole fractions against.
Definition ct_defs.h:173
double takahashiCorrectionFactor(double Pr, double Tr)
Returns interpolated value of (DP)_R obtained from the data in Table 2 of the Takahashi 1975 paper,...
double neufeldCollisionIntegral(double T_star)
Returns the value of the Neufeld collision integral for a given dimensionless temperature.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...