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