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