21 #define COLL_INT_POLY_DEGREE 8
23 GasTransport::GasTransport() :
28 void GasTransport::update_T()
30 if (m_thermo->nSpecies() != m_nsp) {
32 init(m_thermo, m_mode, m_log_level);
35 double T = m_thermo->temperature();
43 m_sqrt_t = sqrt(m_temp);
44 m_t14 = sqrt(m_sqrt_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;
60 double GasTransport::viscosity()
75 multiply(m_phi, m_molefracs.data(), m_spwork.data());
77 for (
size_t k = 0; k < m_nsp; k++) {
78 vismix += m_molefracs[k] * m_visc[k]/m_spwork[k];
84 void GasTransport::updateViscosity_T()
87 updateSpeciesViscosities();
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];
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);
105 void GasTransport::updateSpeciesViscosities()
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]);
114 for (
size_t k = 0; k < m_nsp; k++) {
116 m_sqvisc[k] = m_t14 *
dot5(m_polytempvec, m_visccoeffs[k]);
117 m_visc[k] = (m_sqvisc[k] * m_sqvisc[k]);
123 void GasTransport::updateDiff_T()
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);
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,
141 m_bdiff(j,i) = m_bdiff(i,j);
149 void GasTransport::getBinaryDiffCoeffs(
const size_t ld,
double*
const d)
157 throw CanteraError(
"GasTransport::getBinaryDiffCoeffs",
"ld is too small");
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);
167 void GasTransport::getMixDiffCoeffs(
double*
const d)
177 double mmw = m_thermo->meanMolecularWeight();
178 double p = m_thermo->pressure();
180 d[0] = m_bdiff(0,0) / p;
182 for (
size_t k = 0; k < m_nsp; k++) {
184 for (
size_t j = 0; j < m_nsp; j++) {
186 sum2 += m_molefracs[j] / m_bdiff(j,k);
190 d[k] = m_bdiff(k,k) / p;
192 d[k] = (mmw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
198 void GasTransport::getMixDiffCoeffsMole(
double*
const d)
208 double p = m_thermo->pressure();
210 d[0] = m_bdiff(0,0) / p;
212 for (
size_t k = 0; k < m_nsp; k++) {
214 for (
size_t j = 0; j < m_nsp; j++) {
216 sum2 += m_molefracs[j] / m_bdiff(j,k);
220 d[k] = m_bdiff(k,k) / p;
222 d[k] = (1 - m_molefracs[k]) / (p * sum2);
228 void GasTransport::getMixDiffCoeffsMass(
double*
const d)
238 double mmw = m_thermo->meanMolecularWeight();
239 double p = m_thermo->pressure();
242 d[0] = m_bdiff(0,0) / p;
244 for (
size_t k=0; k<m_nsp; k++) {
247 for (
size_t i=0; i<m_nsp; i++) {
251 sum1 += m_molefracs[i] / m_bdiff(k,i);
252 sum2 += m_molefracs[i] * m_mw[i] / m_bdiff(k,i);
255 sum2 *= p * m_molefracs[k] / (mmw - m_mw[k]*m_molefracs[k]);
256 d[k] = 1.0 / (sum1 + sum2);
261 void GasTransport::init(
ThermoPhase* thermo,
int mode,
int log_level)
266 m_log_level = log_level;
269 setupCollisionParameters();
270 setupCollisionIntegral();
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);
280 m_mw = m_thermo->molecularWeights();
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]);
293 void GasTransport::setupCollisionParameters()
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);
308 m_w_ac.resize(m_nsp);
309 m_disp.resize(m_nsp, 0.0);
310 m_quad_polar.resize(m_nsp, 0.0);
312 const vector<double>& mw = m_thermo->molecularWeights();
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);
320 double f_eps, f_sigma;
322 for (
size_t i = 0; i < m_nsp; i++) {
323 for (
size_t j = i; j < m_nsp; j++) {
325 m_reducedMass(i,j) = mw[i] * mw[j] / (
Avogadro * (mw[i] + mw[j]));
328 m_diam(i,j) = 0.5*(m_sigma[i] + m_sigma[j]);
331 m_epsilon(i,j) = sqrt(m_eps[i]*m_eps[j]);
334 m_dipole(i,j) = sqrt(m_dipole(i,i)*m_dipole(j,j));
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;
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);
354 void GasTransport::setupCollisionIntegral()
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++) {
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));
367 if (m_mode == CK_Mode) {
373 debuglog(
"*** collision_integrals ***\n", m_log_level);
375 integrals.
init(tstar_min, tstar_max, m_log_level);
376 fitCollisionIntegrals(integrals);
377 debuglog(
"*** end of collision_integrals ***\n", m_log_level);
379 debuglog(
"*** property fits ***\n", m_log_level);
380 fitProperties(integrals);
381 debuglog(
"*** end of property fits ***\n", m_log_level);
384 void GasTransport::getTransportData()
386 for (
size_t k = 0; k < m_thermo->nSpecies(); k++) {
387 shared_ptr<Species> s = m_thermo->species(k);
392 "Missing gas-phase transport data for species '{}'.", s->name);
397 }
else if (sptran->
geometry ==
"linear") {
399 }
else if (sptran->
geometry ==
"nonlinear") {
405 m_dipole(k,k) = sptran->
dipole;
406 m_polar[k] = (sptran->
dipole > 0);
409 if (s->input.hasKey(
"critical-parameters") &&
410 s->input[
"critical-parameters"].hasKey(
"acentric-factor"))
412 m_w_ac[k] = s->input[
"critical-parameters"][
"acentric-factor"].asDouble();
421 void GasTransport::makePolarCorrections(
size_t i,
size_t j,
422 double& f_eps,
double& f_sigma)
425 if (m_polar[i] == m_polar[j]) {
433 size_t kp = (m_polar[i] ? i : j);
434 size_t knp = (i == kp ? j : i);
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);
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");
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++) {
466 double dstar = (m_mode != CK_Mode) ? m_delta(i,j) : 0.0;
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);
488 m_poly[i][j] =
static_cast<int>((dptr - fitlist.begin()));
489 m_star_poly_uses_actualT[i][j] = 0;
491 m_poly[j][i] = m_poly[i][j];
492 m_star_poly_uses_actualT[j][i] = m_star_poly_uses_actualT[i][j];
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);
506 m_visccoeffs.clear();
507 m_condcoeffs.clear();
510 for (
size_t n = 0; n < np; n++) {
511 double t = m_thermo->minTemp() + dt*n;
516 vector<double> c(degree + 1), c2(degree + 1);
519 if (m_log_level && m_log_level < 2) {
520 writelog(
"*** polynomial coefficients not printed (log_level < 2) ***\n");
522 double visc, err, relerr,
523 mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
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");
530 writelogf(
"viscosity/sqrt(T) fit to polynomial of degree "
531 "%d in log(T)", degree);
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];
541 double fz_298 = 1.0 + pow(
Pi, 1.5) / sqrt(tstar) * (0.5 + 1.0 / tstar) +
542 (0.25 *
Pi *
Pi + 2) / tstar;
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];
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));
556 double diffcoeff = 3.0/16.0 * sqrt(2.0 *
Pi/m_reducedMass(k,k)) *
558 (
Pi * m_sigma[k] * m_sigma[k] * om11);
562 (om22 *
Pi * m_sigma[k]*m_sigma[k]);
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);
578 if (m_mode == CK_Mode) {
579 spvisc[n] = log(visc);
580 spcond[n] = log(cond);
590 spvisc[n] = sqrt(visc/sqrt_T);
595 spcond[n] = cond/sqrt_T;
596 w[n] = 1.0/(spvisc[n]*spvisc[n]);
597 w2[n] = 1.0/(spcond[n]*spcond[n]);
600 polyfit(np, degree, tlog.data(), spvisc.data(), w.data(), c.data());
601 polyfit(np, degree, tlog.data(), spcond.data(), w2.data(), c2.data());
604 for (
size_t n = 0; n < np; n++) {
606 if (m_mode == CK_Mode) {
607 val = exp(spvisc[n]);
608 fit = exp(
poly3(tlog[n], c.data()));
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);
616 mxerr = std::max(mxerr, fabs(err));
617 mxrelerr = std::max(mxrelerr, fabs(relerr));
621 for (
size_t n = 0; n < np; n++) {
623 if (m_mode == CK_Mode) {
624 val = exp(spcond[n]);
625 fit = exp(
poly3(tlog[n], c2.data()));
627 double sqrt_T = exp(0.5*tlog[n]);
628 val = sqrt_T * spcond[n];
629 fit = sqrt_T *
poly4(tlog[n], c2.data());
633 mxerr_cond = std::max(mxerr_cond, fabs(err));
634 mxrelerr_cond = std::max(mxrelerr_cond, fabs(relerr));
636 m_visccoeffs.push_back(c);
637 m_condcoeffs.push_back(c2);
639 if (m_log_level >= 2) {
643 m_thermo->setTemperature(T_save);
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)");
653 "polynomial of degree %d in log(T)", degree);
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");
661 writelogf(
"Maximum conductivity absolute error: %12.6g\n", mxerr_cond);
662 writelogf(
"Maximum conductivity relative error: %12.6g\n", mxrelerr_cond);
665 writelogf(
"\nbinary diffusion coefficients:\n");
666 if (m_mode == CK_Mode) {
667 writelog(
"log(D) fit to cubic polynomial in log(T)");
669 writelogf(
"D/T**(3/2) fit to polynomial of degree %d in log(T)",degree);
673 fitDiffCoeffs(integrals);
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);
686 for (
size_t n = 0; n < np; n++) {
687 double t = m_thermo->minTemp() + dt*n;
692 vector<double> c(degree + 1), c2(degree + 1);
694 mxerr = 0.0, mxrelerr = 0.0;
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);
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);
712 getBinDiffCorrection(t, integrals, k, j, 1.0, 1.0, fkj, fjk);
714 if (m_mode == CK_Mode) {
715 diff[n] = log(diffcoeff);
718 diff[n] = diffcoeff/pow(t, 1.5);
719 w[n] = 1.0/(diff[n]*diff[n]);
722 polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data());
724 for (
size_t n = 0; n < np; n++) {
726 if (m_mode == CK_Mode) {
728 fit = exp(
poly3(tlog[n], c.data()));
730 double t = exp(tlog[n]);
731 double pre = pow(t, 1.5);
733 fit = pre *
poly4(tlog[n], c.data());
737 mxerr = std::max(mxerr, fabs(err));
738 mxrelerr = std::max(mxrelerr, fabs(relerr));
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");
748 writelogf(
"Maximum binary diffusion coefficient absolute error:"
750 writelogf(
"Maximum binary diffusion coefficient relative error:"
756 size_t k,
size_t j,
double xk,
double xj,
double& fkj,
double& fjk)
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));
779 double cnst = sigratio * sqrt(2.0*w2/wsum) * 2.0 * w1*w1/(wsum * w2);
780 double p1 = cnst * om22_1 / om11_12;
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);
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);
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;
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);
807 void GasTransport::getViscosityPolynomial(
size_t i,
double* coeffs)
const
809 for (
int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
810 coeffs[k] = m_visccoeffs[i][k];
814 void GasTransport::getConductivityPolynomial(
size_t i,
double* coeffs)
const
816 for (
int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
817 coeffs[k] = m_condcoeffs[i][k];
821 void GasTransport::getBinDiffusivityPolynomial(
size_t i,
size_t j,
double* coeffs)
const
823 size_t mi = (j >= i? i : j);
824 size_t mj = (j >= i? j : i);
826 for (
size_t ii = 0; ii < mi; ii++) {
831 for (
int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
832 coeffs[k] = m_diffcoeffs[ic][k];
836 void GasTransport::getCollisionIntegralPolynomial(
size_t i,
size_t j,
837 double* astar_coeffs,
838 double* bstar_coeffs,
839 double* cstar_coeffs)
const
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];
848 void GasTransport::setViscosityPolynomial(
size_t i,
double* coeffs)
850 for (
int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
851 m_visccoeffs[i][k] = coeffs[k];
857 m_bindiff_ok =
false;
861 void GasTransport::setConductivityPolynomial(
size_t i,
double* coeffs)
863 for (
int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
864 m_condcoeffs[i][k] = coeffs[k];
870 m_bindiff_ok =
false;
874 void GasTransport::setBinDiffusivityPolynomial(
size_t i,
size_t j,
double* coeffs)
876 size_t mi = (j >= i? i : j);
877 size_t mj = (j >= i? j : i);
879 for (
size_t ii = 0; ii < mi; ii++) {
884 for (
int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
885 m_diffcoeffs[ic][k] = coeffs[k];
891 m_bindiff_ok =
false;
895 void GasTransport::setCollisionIntegralPolynomial(
size_t i,
size_t j,
896 double* astar_coeffs,
897 double* bstar_coeffs,
898 double* cstar_coeffs,
bool actualT)
901 vector<double> ca(degree+1), cb(degree+1), cc(degree+1);
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];
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];
915 m_star_poly_uses_actualT[i][j] = 1;
916 m_star_poly_uses_actualT[j][i] = m_star_poly_uses_actualT[i][j];
922 m_bindiff_ok =
false;
#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.
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.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
double dot5(const V &x, const V &y)
Templated Inner product of two vectors of length 5.
R poly3(D x, R *c)
Templated evaluation of a polynomial of order 3.
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
double dot4(const V &x, const V &y)
Templated Inner product of two vectors of length 4.
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.
const double Boltzmann
Boltzmann constant [J/K].
const double Avogadro
Avogadro's Number [number/kmol].
const double epsilon_0
Permittivity of free space [F/m].
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
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...