Cantera  3.1.0a1
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"
12 #include "cantera/thermo/Species.h"
13 #include "cantera/base/utilities.h"
14 #include "cantera/base/global.h"
15 
16 namespace 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 
23 GasTransport::GasTransport() :
24  m_polytempvec(5)
25 {
26 }
27 
28 void 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 
60 double 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 
84 void 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 
105 void 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 
123 void 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 
149 void 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 
167 void 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 
198 void 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 
228 void 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 
261 void 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 
293 void 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 
354 void 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 
384 void 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 
421 void 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 
446 void 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 
497 void 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 
676 void 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 
755 void 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 
807 void 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 
814 void 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 
821 void 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 
836 void 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 
848 void 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 
861 void 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 
874 void 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 
895 void 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.
Definition: ctexceptions.h:66
Transport data for a single gas-phase species which can be used in mixture-averaged or multicomponent...
Definition: TransportData.h:46
double polarizability
The polarizability of the molecule [m^3]. Default 0.0.
Definition: TransportData.h:88
double diameter
The Lennard-Jones collision diameter [m].
Definition: TransportData.h:79
double acentric_factor
Pitzer's acentric factor [dimensionless]. Default 0.0.
Definition: TransportData.h:96
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:93
double dispersion_coefficient
dispersion normalized by e^2. [m^5] Default 0.0.
Definition: TransportData.h:99
double dipole
The permanent dipole moment of the molecule [Coulomb-m]. Default 0.0.
Definition: TransportData.h:85
double well_depth
The Lennard-Jones well depth [J].
Definition: TransportData.h:82
string geometry
A string specifying the molecular geometry.
Definition: TransportData.h:76
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.
Definition: ThermoPhase.h:390
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)
Definition: stringUtils.cpp:33
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 poly3(D x, R *c)
Templated evaluation of a polynomial of order 3.
Definition: utilities.h:165
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition: utilities.h:153
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...