Cantera  3.1.0a1
Loading...
Searching...
No Matches
GasTransport.cpp
Go to the documentation of this file.
1//! @file GasTransport.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
7#include "MMCollisionInt.h"
14#include "cantera/base/global.h"
15
16namespace Cantera
17{
18
19//! polynomial degree used for fitting collision integrals
20//! except in CK mode, where the degree is 6.
21#define COLL_INT_POLY_DEGREE 8
22
23GasTransport::GasTransport() :
24 m_polytempvec(5)
25{
26}
27
28void GasTransport::update_T()
29{
30 if (m_thermo->nSpecies() != m_nsp) {
31 // Rebuild data structures if number of species has changed
32 init(m_thermo, m_mode, m_log_level);
33 }
34
35 double T = m_thermo->temperature();
36 if (T == m_temp) {
37 return;
38 }
39
40 m_temp = T;
41 m_kbt = Boltzmann * m_temp;
42 m_logt = log(m_temp);
43 m_sqrt_t = sqrt(m_temp);
44 m_t14 = sqrt(m_sqrt_t);
45
46 // compute powers of log(T)
47 m_polytempvec[0] = 1.0;
48 m_polytempvec[1] = m_logt;
49 m_polytempvec[2] = m_logt*m_logt;
50 m_polytempvec[3] = m_logt*m_logt*m_logt;
51 m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
52
53 // temperature has changed, so polynomial fits will need to be redone
54 m_visc_ok = false;
55 m_spvisc_ok = false;
56 m_viscwt_ok = false;
57 m_bindiff_ok = false;
58}
59
60double GasTransport::viscosity()
61{
62 update_T();
63 update_C();
64
65 if (m_visc_ok) {
66 return m_viscmix;
67 }
68
69 double vismix = 0.0;
70 // update m_visc and m_phi if necessary
71 if (!m_viscwt_ok) {
72 updateViscosity_T();
73 }
74
75 multiply(m_phi, m_molefracs.data(), m_spwork.data());
76
77 for (size_t k = 0; k < m_nsp; k++) {
78 vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
79 }
80 m_viscmix = vismix;
81 return vismix;
82}
83
84void GasTransport::updateViscosity_T()
85{
86 if (!m_spvisc_ok) {
87 updateSpeciesViscosities();
88 }
89
90 // see Eq. (9-5.14) of Poling et al. (2001)
91 for (size_t j = 0; j < m_nsp; j++) {
92 for (size_t k = j; k < m_nsp; k++) {
93 double vratiokj = m_visc[k]/m_visc[j];
94 double wratiojk = m_mw[j]/m_mw[k];
95
96 // Note that m_wratjk(k,j) holds the square root of m_wratjk(j,k)!
97 double factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
98 m_phi(k,j) = factor1*factor1 / (sqrt(8.0) * m_wratkj1(j,k));
99 m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
100 }
101 }
102 m_viscwt_ok = true;
103}
104
105void GasTransport::updateSpeciesViscosities()
106{
107 update_T();
108 if (m_mode == CK_Mode) {
109 for (size_t k = 0; k < m_nsp; k++) {
110 m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
111 m_sqvisc[k] = sqrt(m_visc[k]);
112 }
113 } else {
114 for (size_t k = 0; k < m_nsp; k++) {
115 // the polynomial fit is done for sqrt(visc/sqrt(T))
116 m_sqvisc[k] = m_t14 * dot5(m_polytempvec, m_visccoeffs[k]);
117 m_visc[k] = (m_sqvisc[k] * m_sqvisc[k]);
118 }
119 }
120 m_spvisc_ok = true;
121}
122
123void GasTransport::updateDiff_T()
124{
125 update_T();
126 // evaluate binary diffusion coefficients at unit pressure
127 size_t ic = 0;
128 if (m_mode == CK_Mode) {
129 for (size_t i = 0; i < m_nsp; i++) {
130 for (size_t j = i; j < m_nsp; j++) {
131 m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic]));
132 m_bdiff(j,i) = m_bdiff(i,j);
133 ic++;
134 }
135 }
136 } else {
137 for (size_t i = 0; i < m_nsp; i++) {
138 for (size_t j = i; j < m_nsp; j++) {
139 m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec,
140 m_diffcoeffs[ic]);
141 m_bdiff(j,i) = m_bdiff(i,j);
142 ic++;
143 }
144 }
145 }
146 m_bindiff_ok = true;
147}
148
149void GasTransport::getBinaryDiffCoeffs(const size_t ld, double* const d)
150{
151 update_T();
152 // if necessary, evaluate the binary diffusion coefficients from the polynomial fits
153 if (!m_bindiff_ok) {
154 updateDiff_T();
155 }
156 if (ld < m_nsp) {
157 throw CanteraError("GasTransport::getBinaryDiffCoeffs", "ld is too small");
158 }
159 double rp = 1.0/m_thermo->pressure();
160 for (size_t i = 0; i < m_nsp; i++) {
161 for (size_t j = 0; j < m_nsp; j++) {
162 d[ld*j + i] = rp * m_bdiff(i,j);
163 }
164 }
165}
166
167void GasTransport::getMixDiffCoeffs(double* const d)
168{
169 update_T();
170 update_C();
171
172 // update the binary diffusion coefficients if necessary
173 if (!m_bindiff_ok) {
174 updateDiff_T();
175 }
176
177 double mmw = m_thermo->meanMolecularWeight();
178 double p = m_thermo->pressure();
179 if (m_nsp == 1) {
180 d[0] = m_bdiff(0,0) / p;
181 } else {
182 for (size_t k = 0; k < m_nsp; k++) {
183 double sum2 = 0.0;
184 for (size_t j = 0; j < m_nsp; j++) {
185 if (j != k) {
186 sum2 += m_molefracs[j] / m_bdiff(j,k);
187 }
188 }
189 if (sum2 <= 0.0) {
190 d[k] = m_bdiff(k,k) / p;
191 } else {
192 d[k] = (mmw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
193 }
194 }
195 }
196}
197
198void GasTransport::getMixDiffCoeffsMole(double* const d)
199{
200 update_T();
201 update_C();
202
203 // update the binary diffusion coefficients if necessary
204 if (!m_bindiff_ok) {
205 updateDiff_T();
206 }
207
208 double p = m_thermo->pressure();
209 if (m_nsp == 1) {
210 d[0] = m_bdiff(0,0) / p;
211 } else {
212 for (size_t k = 0; k < m_nsp; k++) {
213 double sum2 = 0.0;
214 for (size_t j = 0; j < m_nsp; j++) {
215 if (j != k) {
216 sum2 += m_molefracs[j] / m_bdiff(j,k);
217 }
218 }
219 if (sum2 <= 0.0) {
220 d[k] = m_bdiff(k,k) / p;
221 } else {
222 d[k] = (1 - m_molefracs[k]) / (p * sum2);
223 }
224 }
225 }
226}
227
228void GasTransport::getMixDiffCoeffsMass(double* const d)
229{
230 update_T();
231 update_C();
232
233 // update the binary diffusion coefficients if necessary
234 if (!m_bindiff_ok) {
235 updateDiff_T();
236 }
237
238 double mmw = m_thermo->meanMolecularWeight();
239 double p = m_thermo->pressure();
240
241 if (m_nsp == 1) {
242 d[0] = m_bdiff(0,0) / p;
243 } else {
244 for (size_t k=0; k<m_nsp; k++) {
245 double sum1 = 0.0;
246 double sum2 = 0.0;
247 for (size_t i=0; i<m_nsp; i++) {
248 if (i==k) {
249 continue;
250 }
251 sum1 += m_molefracs[i] / m_bdiff(k,i);
252 sum2 += m_molefracs[i] * m_mw[i] / m_bdiff(k,i);
253 }
254 sum1 *= p;
255 sum2 *= p * m_molefracs[k] / (mmw - m_mw[k]*m_molefracs[k]);
256 d[k] = 1.0 / (sum1 + sum2);
257 }
258 }
259}
260
261void GasTransport::init(ThermoPhase* thermo, int mode, int log_level)
262{
263 m_thermo = thermo;
264 m_nsp = m_thermo->nSpecies();
265 m_mode = mode;
266 m_log_level = log_level;
267
268 // set up Monchick and Mason collision integrals
269 setupCollisionParameters();
270 setupCollisionIntegral();
271
272 m_molefracs.resize(m_nsp);
273 m_spwork.resize(m_nsp);
274 m_visc.resize(m_nsp);
275 m_sqvisc.resize(m_nsp);
276 m_phi.resize(m_nsp, m_nsp, 0.0);
277 m_bdiff.resize(m_nsp, m_nsp);
278
279 // make a local copy of the molecular weights
280 m_mw = m_thermo->molecularWeights();
281
282 m_wratjk.resize(m_nsp, m_nsp, 0.0);
283 m_wratkj1.resize(m_nsp, m_nsp, 0.0);
284 for (size_t j = 0; j < m_nsp; j++) {
285 for (size_t k = j; k < m_nsp; k++) {
286 m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
287 m_wratjk(k,j) = sqrt(m_wratjk(j,k));
288 m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
289 }
290 }
291}
292
293void GasTransport::setupCollisionParameters()
294{
295 m_epsilon.resize(m_nsp, m_nsp, 0.0);
296 m_delta.resize(m_nsp, m_nsp, 0.0);
297 m_reducedMass.resize(m_nsp, m_nsp, 0.0);
298 m_dipole.resize(m_nsp, m_nsp, 0.0);
299 m_diam.resize(m_nsp, m_nsp, 0.0);
300 m_crot.resize(m_nsp);
301 m_zrot.resize(m_nsp);
302 m_polar.resize(m_nsp, false);
303 m_alpha.resize(m_nsp, 0.0);
304 m_poly.resize(m_nsp);
305 m_star_poly_uses_actualT.resize(m_nsp);
306 m_sigma.resize(m_nsp);
307 m_eps.resize(m_nsp);
308 m_w_ac.resize(m_nsp);
309 m_disp.resize(m_nsp, 0.0);
310 m_quad_polar.resize(m_nsp, 0.0);
311
312 const vector<double>& mw = m_thermo->molecularWeights();
313 getTransportData();
314
315 for (size_t i = 0; i < m_nsp; i++) {
316 m_poly[i].resize(m_nsp);
317 m_star_poly_uses_actualT[i].resize(m_nsp);
318 }
319
320 double f_eps, f_sigma;
321
322 for (size_t i = 0; i < m_nsp; i++) {
323 for (size_t j = i; j < m_nsp; j++) {
324 // the reduced mass
325 m_reducedMass(i,j) = mw[i] * mw[j] / (Avogadro * (mw[i] + mw[j]));
326
327 // hard-sphere diameter for (i,j) collisions
328 m_diam(i,j) = 0.5*(m_sigma[i] + m_sigma[j]);
329
330 // the effective well depth for (i,j) collisions
331 m_epsilon(i,j) = sqrt(m_eps[i]*m_eps[j]);
332
333 // the effective dipole moment for (i,j) collisions
334 m_dipole(i,j) = sqrt(m_dipole(i,i)*m_dipole(j,j));
335
336 // reduced dipole moment delta* (nondimensional)
337 double d = m_diam(i,j);
338 m_delta(i,j) = 0.5 * m_dipole(i,j)*m_dipole(i,j)
339 / (4 * Pi * epsilon_0 * m_epsilon(i,j) * d * d * d);
340 makePolarCorrections(i, j, f_eps, f_sigma);
341 m_diam(i,j) *= f_sigma;
342 m_epsilon(i,j) *= f_eps;
343
344 // properties are symmetric
345 m_reducedMass(j,i) = m_reducedMass(i,j);
346 m_diam(j,i) = m_diam(i,j);
347 m_epsilon(j,i) = m_epsilon(i,j);
348 m_dipole(j,i) = m_dipole(i,j);
349 m_delta(j,i) = m_delta(i,j);
350 }
351 }
352}
353
354void GasTransport::setupCollisionIntegral()
355{
356 double tstar_min = 1.e8, tstar_max = 0.0;
357 for (size_t i = 0; i < m_nsp; i++) {
358 for (size_t j = i; j < m_nsp; j++) {
359 // The polynomial fits of collision integrals vs. T*
360 // will be done for the T* from tstar_min to tstar_max
361 tstar_min = std::min(tstar_min, Boltzmann * m_thermo->minTemp()/m_epsilon(i,j));
362 tstar_max = std::max(tstar_max, Boltzmann * m_thermo->maxTemp()/m_epsilon(i,j));
363 }
364 }
365 // Chemkin fits the entire T* range in the Monchick and Mason tables,
366 // so modify tstar_min and tstar_max if in Chemkin compatibility mode
367 if (m_mode == CK_Mode) {
368 tstar_min = 0.101;
369 tstar_max = 99.9;
370 }
371
372 // initialize the collision integral calculator for the desired T* range
373 debuglog("*** collision_integrals ***\n", m_log_level);
374 MMCollisionInt integrals;
375 integrals.init(tstar_min, tstar_max, m_log_level);
376 fitCollisionIntegrals(integrals);
377 debuglog("*** end of collision_integrals ***\n", m_log_level);
378 // make polynomial fits
379 debuglog("*** property fits ***\n", m_log_level);
380 fitProperties(integrals);
381 debuglog("*** end of property fits ***\n", m_log_level);
382}
383
384void GasTransport::getTransportData()
385{
386 for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
387 shared_ptr<Species> s = m_thermo->species(k);
388 const GasTransportData* sptran =
389 dynamic_cast<GasTransportData*>(s->transport.get());
390 if (!sptran) {
391 throw CanteraError("GasTransport::getTransportData",
392 "Missing gas-phase transport data for species '{}'.", s->name);
393 }
394
395 if (sptran->geometry == "atom") {
396 m_crot[k] = 0.0;
397 } else if (sptran->geometry == "linear") {
398 m_crot[k] = 1.0;
399 } else if (sptran->geometry == "nonlinear") {
400 m_crot[k] = 1.5;
401 }
402
403 m_sigma[k] = sptran->diameter;
404 m_eps[k] = sptran->well_depth;
405 m_dipole(k,k) = sptran->dipole;
406 m_polar[k] = (sptran->dipole > 0);
407 m_alpha[k] = sptran->polarizability;
408 m_zrot[k] = sptran->rotational_relaxation;
409 if (s->input.hasKey("critical-parameters") &&
410 s->input["critical-parameters"].hasKey("acentric-factor"))
411 {
412 m_w_ac[k] = s->input["critical-parameters"]["acentric-factor"].asDouble();
413 } else {
414 m_w_ac[k] = sptran->acentric_factor;
415 }
416 m_disp[k] = sptran->dispersion_coefficient;
417 m_quad_polar[k] = sptran->quadrupole_polarizability;
418 }
419}
420
421void GasTransport::makePolarCorrections(size_t i, size_t j,
422 double& f_eps, double& f_sigma)
423{
424 // no correction if both are nonpolar, or both are polar
425 if (m_polar[i] == m_polar[j]) {
426 f_eps = 1.0;
427 f_sigma = 1.0;
428 return;
429 }
430
431 // corrections to the effective diameter and well depth
432 // if one is polar and one is non-polar
433 size_t kp = (m_polar[i] ? i : j); // the polar one
434 size_t knp = (i == kp ? j : i); // the nonpolar one
435 double d3np, d3p, alpha_star, mu_p_star, xi;
436 d3np = pow(m_sigma[knp],3);
437 d3p = pow(m_sigma[kp],3);
438 alpha_star = m_alpha[knp]/d3np;
439 mu_p_star = m_dipole(kp,kp)/sqrt(4 * Pi * epsilon_0 * d3p * m_eps[kp]);
440 xi = 1.0 + 0.25 * alpha_star * mu_p_star * mu_p_star *
441 sqrt(m_eps[kp]/m_eps[knp]);
442 f_sigma = pow(xi, -1.0/6.0);
443 f_eps = xi*xi;
444}
445
446void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals)
447{
448 // Chemkin fits to sixth order polynomials
449 int degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
450 if (m_log_level) {
451 writelog("tstar_fits\n"
452 "fits to A*, B*, and C* vs. log(T*).\n"
453 "These are done only for the required dstar(j,k) values.\n\n");
454 if (m_log_level < 3) {
455 writelog("*** polynomial coefficients not printed (log_level < 3) ***\n");
456 }
457 }
458 vector<double> fitlist;
459 m_omega22_poly.clear();
460 m_astar_poly.clear();
461 m_bstar_poly.clear();
462 m_cstar_poly.clear();
463 for (size_t i = 0; i < m_nsp; i++) {
464 for (size_t j = i; j < m_nsp; j++) {
465 // Chemkin fits only delta* = 0
466 double dstar = (m_mode != CK_Mode) ? m_delta(i,j) : 0.0;
467
468 // if a fit has already been generated for delta* = m_delta(i,j),
469 // then use it. Otherwise, make a new fit, and add m_delta(i,j) to
470 // the list of delta* values for which fits have been done.
471
472 // 'find' returns a pointer to end() if not found
473 auto dptr = find(fitlist.begin(), fitlist.end(), dstar);
474 if (dptr == fitlist.end()) {
475 vector<double> ca(degree+1), cb(degree+1), cc(degree+1);
476 vector<double> co22(degree+1);
477 integrals.fit(degree, dstar, ca.data(), cb.data(), cc.data());
478 integrals.fit_omega22(degree, dstar, co22.data());
479 m_omega22_poly.push_back(co22);
480 m_astar_poly.push_back(ca);
481 m_bstar_poly.push_back(cb);
482 m_cstar_poly.push_back(cc);
483 m_poly[i][j] = static_cast<int>(m_astar_poly.size()) - 1;
484 m_star_poly_uses_actualT[i][j] = 0;
485 fitlist.push_back(dstar);
486 } else {
487 // delta* found in fitlist, so just point to this polynomial
488 m_poly[i][j] = static_cast<int>((dptr - fitlist.begin()));
489 m_star_poly_uses_actualT[i][j] = 0;
490 }
491 m_poly[j][i] = m_poly[i][j];
492 m_star_poly_uses_actualT[j][i] = m_star_poly_uses_actualT[i][j];
493 }
494 }
495}
496
497void GasTransport::fitProperties(MMCollisionInt& integrals)
498{
499 // number of points to use in generating fit data
500 const size_t np = 50;
501 int degree = (m_mode == CK_Mode ? 3 : 4);
502 double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1);
503 vector<double> tlog(np), spvisc(np), spcond(np);
504 vector<double> w(np), w2(np);
505
506 m_visccoeffs.clear();
507 m_condcoeffs.clear();
508
509 // generate array of log(t) values
510 for (size_t n = 0; n < np; n++) {
511 double t = m_thermo->minTemp() + dt*n;
512 tlog[n] = log(t);
513 }
514
515 // vector of polynomial coefficients
516 vector<double> c(degree + 1), c2(degree + 1);
517
518 // fit the pure-species viscosity and thermal conductivity for each species
519 if (m_log_level && m_log_level < 2) {
520 writelog("*** polynomial coefficients not printed (log_level < 2) ***\n");
521 }
522 double visc, err, relerr,
523 mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
524
525 if (m_log_level) {
526 writelog("Polynomial fits for viscosity:\n");
527 if (m_mode == CK_Mode) {
528 writelog("log(viscosity) fit to cubic polynomial in log(T)\n");
529 } else {
530 writelogf("viscosity/sqrt(T) fit to polynomial of degree "
531 "%d in log(T)", degree);
532 }
533 }
534
535 double T_save = m_thermo->temperature();
536 const vector<double>& mw = m_thermo->molecularWeights();
537 for (size_t k = 0; k < m_nsp; k++) {
538 double tstar = Boltzmann * 298.0 / m_eps[k];
539 // Scaling factor for temperature dependence of z_rot. [Kee2003] Eq.
540 // 12.112 or [Kee2017] Eq. 11.115
541 double fz_298 = 1.0 + pow(Pi, 1.5) / sqrt(tstar) * (0.5 + 1.0 / tstar) +
542 (0.25 * Pi * Pi + 2) / tstar;
543
544 for (size_t n = 0; n < np; n++) {
545 double t = m_thermo->minTemp() + dt*n;
546 m_thermo->setTemperature(t);
547 vector<double> cp_R_all(m_thermo->nSpecies());
548 m_thermo->getCp_R_ref(&cp_R_all[0]);
549 double cp_R = cp_R_all[k];
550 tstar = Boltzmann * t / m_eps[k];
551 double sqrt_T = sqrt(t);
552 double om22 = integrals.omega22(tstar, m_delta(k,k));
553 double om11 = integrals.omega11(tstar, m_delta(k,k));
554
555 // self-diffusion coefficient, without polar corrections
556 double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,k)) *
557 pow((Boltzmann * t), 1.5)/
558 (Pi * m_sigma[k] * m_sigma[k] * om11);
559
560 // viscosity
561 visc = 5.0/16.0 * sqrt(Pi * mw[k] * Boltzmann * t / Avogadro) /
562 (om22 * Pi * m_sigma[k]*m_sigma[k]);
563
564 // thermal conductivity
565 double f_int = mw[k]/(GasConstant * t) * diffcoeff/visc;
566 double cv_rot = m_crot[k];
567 double A_factor = 2.5 - f_int;
568 double fz_tstar = 1.0 + pow(Pi, 1.5) / sqrt(tstar) * (0.5 + 1.0 / tstar) +
569 (0.25 * Pi * Pi + 2) / tstar;
570 double B_factor = m_zrot[k] * fz_298 / fz_tstar + 2.0/Pi * (5.0/3.0 * cv_rot + f_int);
571 double c1 = 2.0/Pi * A_factor/B_factor;
572 double cv_int = cp_R - 2.5 - cv_rot;
573 double f_rot = f_int * (1.0 + c1);
574 double f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5);
575 double cond = (visc/mw[k])*GasConstant*(f_trans * 1.5
576 + f_rot * cv_rot + f_int * cv_int);
577
578 if (m_mode == CK_Mode) {
579 spvisc[n] = log(visc);
580 spcond[n] = log(cond);
581 w[n] = -1.0;
582 w2[n] = -1.0;
583 } else {
584 // the viscosity should be proportional approximately to
585 // sqrt(T); therefore, visc/sqrt(T) should have only a weak
586 // temperature dependence. And since the mixture rule requires
587 // the square root of the pure-species viscosity, fit the square
588 // root of (visc/sqrt(T)) to avoid having to compute square
589 // roots in the mixture rule.
590 spvisc[n] = sqrt(visc/sqrt_T);
591
592 // the pure-species conductivity scales approximately with
593 // sqrt(T). Unlike the viscosity, there is no reason here to fit
594 // the square root, since a different mixture rule is used.
595 spcond[n] = cond/sqrt_T;
596 w[n] = 1.0/(spvisc[n]*spvisc[n]);
597 w2[n] = 1.0/(spcond[n]*spcond[n]);
598 }
599 }
600 polyfit(np, degree, tlog.data(), spvisc.data(), w.data(), c.data());
601 polyfit(np, degree, tlog.data(), spcond.data(), w2.data(), c2.data());
602
603 // evaluate max fit errors for viscosity
604 for (size_t n = 0; n < np; n++) {
605 double val, fit;
606 if (m_mode == CK_Mode) {
607 val = exp(spvisc[n]);
608 fit = exp(poly3(tlog[n], c.data()));
609 } else {
610 double sqrt_T = exp(0.5*tlog[n]);
611 val = sqrt_T * pow(spvisc[n],2);
612 fit = sqrt_T * pow(poly4(tlog[n], c.data()),2);
613 }
614 err = fit - val;
615 relerr = err/val;
616 mxerr = std::max(mxerr, fabs(err));
617 mxrelerr = std::max(mxrelerr, fabs(relerr));
618 }
619
620 // evaluate max fit errors for conductivity
621 for (size_t n = 0; n < np; n++) {
622 double val, fit;
623 if (m_mode == CK_Mode) {
624 val = exp(spcond[n]);
625 fit = exp(poly3(tlog[n], c2.data()));
626 } else {
627 double sqrt_T = exp(0.5*tlog[n]);
628 val = sqrt_T * spcond[n];
629 fit = sqrt_T * poly4(tlog[n], c2.data());
630 }
631 err = fit - val;
632 relerr = err/val;
633 mxerr_cond = std::max(mxerr_cond, fabs(err));
634 mxrelerr_cond = std::max(mxrelerr_cond, fabs(relerr));
635 }
636 m_visccoeffs.push_back(c);
637 m_condcoeffs.push_back(c2);
638
639 if (m_log_level >= 2) {
640 writelog(m_thermo->speciesName(k) + ": [" + vec2str(c) + "]\n");
641 }
642 }
643 m_thermo->setTemperature(T_save);
644
645 if (m_log_level) {
646 writelogf("Maximum viscosity absolute error: %12.6g\n", mxerr);
647 writelogf("Maximum viscosity relative error: %12.6g\n", mxrelerr);
648 writelog("\nPolynomial fits for conductivity:\n");
649 if (m_mode == CK_Mode) {
650 writelog("log(conductivity) fit to cubic polynomial in log(T)");
651 } else {
652 writelogf("conductivity/sqrt(T) fit to "
653 "polynomial of degree %d in log(T)", degree);
654 }
655 if (m_log_level >= 2) {
656 for (size_t k = 0; k < m_nsp; k++) {
657 writelog(m_thermo->speciesName(k) + ": [" +
658 vec2str(m_condcoeffs[k]) + "]\n");
659 }
660 }
661 writelogf("Maximum conductivity absolute error: %12.6g\n", mxerr_cond);
662 writelogf("Maximum conductivity relative error: %12.6g\n", mxrelerr_cond);
663
664 // fit the binary diffusion coefficients for each species pair
665 writelogf("\nbinary diffusion coefficients:\n");
666 if (m_mode == CK_Mode) {
667 writelog("log(D) fit to cubic polynomial in log(T)");
668 } else {
669 writelogf("D/T**(3/2) fit to polynomial of degree %d in log(T)",degree);
670 }
671 }
672
673 fitDiffCoeffs(integrals);
674}
675
676void GasTransport::fitDiffCoeffs(MMCollisionInt& integrals)
677{
678 // number of points to use in generating fit data
679 const size_t np = 50;
680 int degree = (m_mode == CK_Mode ? 3 : 4);
681 double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1);
682 vector<double> tlog(np);
683 vector<double> w(np), w2(np);
684
685 // generate array of log(t) values
686 for (size_t n = 0; n < np; n++) {
687 double t = m_thermo->minTemp() + dt*n;
688 tlog[n] = log(t);
689 }
690
691 // vector of polynomial coefficients
692 vector<double> c(degree + 1), c2(degree + 1);
693 double err, relerr,
694 mxerr = 0.0, mxrelerr = 0.0;
695
696 vector<double> diff(np + 1);
697 m_diffcoeffs.clear();
698 for (size_t k = 0; k < m_nsp; k++) {
699 for (size_t j = k; j < m_nsp; j++) {
700 for (size_t n = 0; n < np; n++) {
701 double t = m_thermo->minTemp() + dt*n;
702 double eps = m_epsilon(j,k);
703 double tstar = Boltzmann * t/eps;
704 double sigma = m_diam(j,k);
705 double om11 = integrals.omega11(tstar, m_delta(j,k));
706 double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j))
707 * pow(Boltzmann * t, 1.5) / (Pi * sigma * sigma * om11);
708
709 // 2nd order correction
710 // NOTE: THIS CORRECTION IS NOT APPLIED
711 double fkj, fjk;
712 getBinDiffCorrection(t, integrals, k, j, 1.0, 1.0, fkj, fjk);
713
714 if (m_mode == CK_Mode) {
715 diff[n] = log(diffcoeff);
716 w[n] = -1.0;
717 } else {
718 diff[n] = diffcoeff/pow(t, 1.5);
719 w[n] = 1.0/(diff[n]*diff[n]);
720 }
721 }
722 polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data());
723
724 for (size_t n = 0; n < np; n++) {
725 double val, fit;
726 if (m_mode == CK_Mode) {
727 val = exp(diff[n]);
728 fit = exp(poly3(tlog[n], c.data()));
729 } else {
730 double t = exp(tlog[n]);
731 double pre = pow(t, 1.5);
732 val = pre * diff[n];
733 fit = pre * poly4(tlog[n], c.data());
734 }
735 err = fit - val;
736 relerr = err/val;
737 mxerr = std::max(mxerr, fabs(err));
738 mxrelerr = std::max(mxrelerr, fabs(relerr));
739 }
740 m_diffcoeffs.push_back(c);
741 if (m_log_level >= 2) {
742 writelog(m_thermo->speciesName(k) + "__" +
743 m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
744 }
745 }
746 }
747 if (m_log_level) {
748 writelogf("Maximum binary diffusion coefficient absolute error:"
749 " %12.6g\n", mxerr);
750 writelogf("Maximum binary diffusion coefficient relative error:"
751 "%12.6g", mxrelerr);
752 }
753}
754
755void GasTransport::getBinDiffCorrection(double t, MMCollisionInt& integrals,
756 size_t k, size_t j, double xk, double xj, double& fkj, double& fjk)
757{
758 double w1 = m_thermo->molecularWeight(k);
759 double w2 = m_thermo->molecularWeight(j);
760 double wsum = w1 + w2;
761 double wmwp = (w1 - w2)/wsum;
762 double sqw12 = sqrt(w1*w2);
763 double sig1 = m_sigma[k];
764 double sig2 = m_sigma[j];
765 double sig12 = 0.5*(m_sigma[k] + m_sigma[j]);
766 double sigratio = sig1*sig1/(sig2*sig2);
767 double sigratio2 = sig1*sig1/(sig12*sig12);
768 double sigratio3 = sig2*sig2/(sig12*sig12);
769 double tstar1 = Boltzmann * t / m_eps[k];
770 double tstar2 = Boltzmann * t / m_eps[j];
771 double tstar12 = Boltzmann * t / sqrt(m_eps[k] * m_eps[j]);
772 double om22_1 = integrals.omega22(tstar1, m_delta(k,k));
773 double om22_2 = integrals.omega22(tstar2, m_delta(j,j));
774 double om11_12 = integrals.omega11(tstar12, m_delta(k,j));
775 double astar_12 = integrals.astar(tstar12, m_delta(k,j));
776 double bstar_12 = integrals.bstar(tstar12, m_delta(k,j));
777 double cstar_12 = integrals.cstar(tstar12, m_delta(k,j));
778
779 double cnst = sigratio * sqrt(2.0*w2/wsum) * 2.0 * w1*w1/(wsum * w2);
780 double p1 = cnst * om22_1 / om11_12;
781
782 cnst = (1.0/sigratio) * sqrt(2.0*w1/wsum) * 2.0*w2*w2/(wsum*w1);
783 double p2 = cnst * om22_2 / om11_12;
784 double p12 = 15.0 * wmwp*wmwp + 8.0*w1*w2*astar_12/(wsum*wsum);
785
786 cnst = (2.0/(w2*wsum))*sqrt(2.0*w2/wsum)*sigratio2;
787 double q1 = cnst*((2.5 - 1.2*bstar_12)*w1*w1 + 3.0*w2*w2
788 + 1.6*w1*w2*astar_12);
789
790 cnst = (2.0/(w1*wsum))*sqrt(2.0*w1/wsum)*sigratio3;
791 double q2 = cnst*((2.5 - 1.2*bstar_12)*w2*w2 + 3.0*w1*w1
792 + 1.6*w1*w2*astar_12);
793 double q12 = wmwp*wmwp*15.0*(2.5 - 1.2*bstar_12)
794 + 4.0*w1*w2*astar_12*(11.0 - 2.4*bstar_12)/(wsum*wsum)
795 + 1.6*wsum*om22_1*om22_2/(om11_12*om11_12*sqw12)
796 * sigratio2 * sigratio3;
797
798 cnst = 6.0*cstar_12 - 5.0;
799 fkj = 1.0 + 0.1*cnst*cnst *
800 (p1*xk*xk + p2*xj*xj + p12*xk*xj)/
801 (q1*xk*xk + q2*xj*xj + q12*xk*xj);
802 fjk = 1.0 + 0.1*cnst*cnst *
803 (p2*xk*xk + p1*xj*xj + p12*xk*xj)/
804 (q2*xk*xk + q1*xj*xj + q12*xk*xj);
805}
806
807void GasTransport::getViscosityPolynomial(size_t i, double* coeffs) const
808{
809 for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
810 coeffs[k] = m_visccoeffs[i][k];
811 }
812}
813
814void GasTransport::getConductivityPolynomial(size_t i, double* coeffs) const
815{
816 for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
817 coeffs[k] = m_condcoeffs[i][k];
818 }
819}
820
821void GasTransport::getBinDiffusivityPolynomial(size_t i, size_t j, double* coeffs) const
822{
823 size_t mi = (j >= i? i : j);
824 size_t mj = (j >= i? j : i);
825 size_t ic = 0;
826 for (size_t ii = 0; ii < mi; ii++) {
827 ic += m_nsp - ii;
828 }
829 ic += mj - mi;
830
831 for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
832 coeffs[k] = m_diffcoeffs[ic][k];
833 }
834}
835
836void GasTransport::getCollisionIntegralPolynomial(size_t i, size_t j,
837 double* astar_coeffs,
838 double* bstar_coeffs,
839 double* cstar_coeffs) const
840{
841 for (int k = 0; k < (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE) + 1; k++) {
842 astar_coeffs[k] = m_astar_poly[m_poly[i][j]][k];
843 bstar_coeffs[k] = m_bstar_poly[m_poly[i][j]][k];
844 cstar_coeffs[k] = m_cstar_poly[m_poly[i][j]][k];
845 }
846}
847
848void GasTransport::setViscosityPolynomial(size_t i, double* coeffs)
849{
850 for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
851 m_visccoeffs[i][k] = coeffs[k];
852 }
853
854 m_visc_ok = false;
855 m_spvisc_ok = false;
856 m_viscwt_ok = false;
857 m_bindiff_ok = false;
858 m_temp = -1;
859}
860
861void GasTransport::setConductivityPolynomial(size_t i, double* coeffs)
862{
863 for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
864 m_condcoeffs[i][k] = coeffs[k];
865 }
866
867 m_visc_ok = false;
868 m_spvisc_ok = false;
869 m_viscwt_ok = false;
870 m_bindiff_ok = false;
871 m_temp = -1;
872}
873
874void GasTransport::setBinDiffusivityPolynomial(size_t i, size_t j, double* coeffs)
875{
876 size_t mi = (j >= i? i : j);
877 size_t mj = (j >= i? j : i);
878 size_t ic = 0;
879 for (size_t ii = 0; ii < mi; ii++) {
880 ic += m_nsp - ii;
881 }
882 ic += mj - mi;
883
884 for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
885 m_diffcoeffs[ic][k] = coeffs[k];
886 }
887
888 m_visc_ok = false;
889 m_spvisc_ok = false;
890 m_viscwt_ok = false;
891 m_bindiff_ok = false;
892 m_temp = -1;
893}
894
895void GasTransport::setCollisionIntegralPolynomial(size_t i, size_t j,
896 double* astar_coeffs,
897 double* bstar_coeffs,
898 double* cstar_coeffs, bool actualT)
899{
900 size_t degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
901 vector<double> ca(degree+1), cb(degree+1), cc(degree+1);
902
903 for (size_t k = 0; k < degree+1; k++) {
904 ca[k] = astar_coeffs[k];
905 cb[k] = bstar_coeffs[k];
906 cc[k] = cstar_coeffs[k];
907 }
908
909 m_astar_poly.push_back(ca);
910 m_bstar_poly.push_back(cb);
911 m_cstar_poly.push_back(cc);
912 m_poly[i][j] = static_cast<int>(m_astar_poly.size()) - 1;
913 m_poly[j][i] = m_poly[i][j];
914 if (actualT) {
915 m_star_poly_uses_actualT[i][j] = 1;
916 m_star_poly_uses_actualT[j][i] = m_star_poly_uses_actualT[i][j];
917 }
918
919 m_visc_ok = false;
920 m_spvisc_ok = false;
921 m_viscwt_ok = false;
922 m_bindiff_ok = false;
923 m_temp = -1;
924}
925
926}
#define COLL_INT_POLY_DEGREE
polynomial degree used for fitting collision integrals except in CK mode, where the degree is 6.
Monchick and Mason collision integrals.
Declaration for class Cantera::Species.
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class for exceptions thrown by Cantera classes.
Transport data for a single gas-phase species which can be used in mixture-averaged or multicomponent...
double polarizability
The polarizability of the molecule [m^3]. Default 0.0.
double diameter
The Lennard-Jones collision diameter [m].
double acentric_factor
Pitzer's acentric factor [dimensionless]. Default 0.0.
double quadrupole_polarizability
quadrupole. Default 0.0.
double rotational_relaxation
The rotational relaxation number (the number of collisions it takes to equilibrate the rotational deg...
double dispersion_coefficient
dispersion normalized by e^2. [m^5] Default 0.0.
double dipole
The permanent dipole moment of the molecule [Coulomb-m]. Default 0.0.
double well_depth
The Lennard-Jones well depth [J].
string geometry
A string specifying the molecular geometry.
Calculation of Collision integrals.
void init(double tsmin, double tsmax, int loglevel=0)
Initialize the object for calculation.
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
Base class for a phase with thermodynamic properties.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
string vec2str(const vector< double > &v, const string &fmt, const string &sep)
Convert a vector to a string (separated by commas)
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition global.h:158
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
double dot5(const V &x, const V &y)
Templated Inner product of two vectors of length 5.
Definition utilities.h:55
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition utilities.h:153
R poly3(D x, R *c)
Templated evaluation of a polynomial of order 3.
Definition utilities.h:165
double dot4(const V &x, const V &y)
Templated Inner product of two vectors of length 4.
Definition utilities.h:40
double polyfit(size_t n, size_t deg, const double *xp, const double *yp, const double *wp, double *pp)
Fits a polynomial function to a set of data points.
Definition polyfit.cpp:14
const double Boltzmann
Boltzmann constant [J/K].
Definition ct_defs.h:84
const double Avogadro
Avogadro's Number [number/kmol].
Definition ct_defs.h:81
const double epsilon_0
Permittivity of free space [F/m].
Definition ct_defs.h:137
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
const double Pi
Pi.
Definition ct_defs.h:68
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...