Cantera  2.0
MultiTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file MultiTransport.cpp
3  * Implementation file for class MultiTransport
4  */
5 /*
6  * Copyright 2001 California Institute of Technology
7  * See file License.txt for licensing information
8  */
9 
11 
14 
16 #include "cantera/base/utilities.h"
17 #include "cantera/base/utilities.h"
18 #include "L_matrix.h"
21 
24 
25 #include <iostream>
26 using namespace std;
27 
28 /**
29  * Mole fractions below MIN_X will be set to MIN_X when computing
30  * transport properties.
31  */
32 #define MIN_X 1.e-20
33 
34 namespace Cantera
35 {
36 
37 ///////////////////// helper functions /////////////////////////
38 
39 /**
40  * @internal
41  *
42  * The Parker temperature correction to the rotational collision
43  * number.
44  *
45  * @param tr Reduced temperature \f$ \epsilon/kT \f$
46  * @param sqtr square root of tr.
47  */
48 inline doublereal Frot(doublereal tr, doublereal sqtr)
49 {
50  const doublereal c1 = 0.5*SqrtPi*Pi;
51  const doublereal c2 = 0.25*Pi*Pi + 2.0;
52  const doublereal c3 = SqrtPi*Pi;
53  return 1.0 + c1*sqtr + c2*tr + c3*sqtr*tr;
54 }
55 
56 
57 /**
58  * This method is used by GMRES to multiply the L matrix by a
59  * vector b. The L matrix has a 3x3 block structure, where each
60  * block is a K x K matrix. The elements of the upper-right and
61  * lower-left blocks are all zero. This method is defined so
62  * that the multiplication only involves the seven non-zero
63  * blocks.
64  */
65 void L_Matrix::mult(const doublereal* b, doublereal* prod) const
66 {
67  integer n = static_cast<int>(nRows())/3;
68  integer n2 = 2*n;
69  integer n3 = 3*n;
70  ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, n, n2, 1.0,
71  DATA_PTR(data()), static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
72  ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, n, n3, 1.0,
73  DATA_PTR(data()) + n, static_cast<int>(nRows()),
74  b, 1, 0.0, prod+n, 1);
75  ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, n, n, 1.0,
76  DATA_PTR(data()) + n*n3 + n2, static_cast<int>(nRows()),
77  b + n, 1, 0.0, prod+n2, 1);
78  for (int i = 0; i < n; i++) {
79  prod[i + n2] += b[i + n2] * value(i + n2, i + n2);
80  }
81 }
82 
83 //////////////////// class MultiTransport methods //////////////
84 
85 MultiTransport::MultiTransport(thermo_t* thermo)
86  : GasTransport(thermo)
87 {
88 }
89 
90 //====================================================================================================================
92 {
94 
95  // copy polynomials and parameters into local storage
96  m_poly = tr.poly;
97  m_astar_poly = tr.astar_poly;
98  m_bstar_poly = tr.bstar_poly;
99  m_cstar_poly = tr.cstar_poly;
100  m_om22_poly = tr.omega22_poly;
101  m_zrot = tr.zrot;
102  m_crot = tr.crot;
103  m_epsilon = tr.epsilon;
104  m_diam = tr.diam;
105  m_eps = tr.eps;
106  m_alpha = tr.alpha;
107  m_dipoleDiag.resize(m_nsp);
108  for (size_t i = 0; i < m_nsp; i++) {
109  m_dipoleDiag[i] = tr.dipole(i,i);
110  }
111 
112  // the L matrix
113  m_Lmatrix.resize(3*m_nsp, 3*m_nsp);
114  m_a.resize(3*m_nsp, 1.0);
115  m_b.resize(3*m_nsp, 0.0);
116  m_aa.resize(m_nsp, m_nsp, 0.0);
117 
118  m_frot_298.resize(m_nsp);
119  m_rotrelax.resize(m_nsp);
120 
121  m_cinternal.resize(m_nsp);
122 
123  m_om22.resize(m_nsp, m_nsp);
124  m_astar.resize(m_nsp, m_nsp);
125  m_bstar.resize(m_nsp, m_nsp);
126  m_cstar.resize(m_nsp, m_nsp);
127 
128  // set flags all false
129  m_abc_ok = false;
130  m_l0000_ok = false;
131  m_lmatrix_soln_ok = false;
132 
133  m_thermal_tlast = 0.0;
134 
135  // use LU decomposition by default
136  m_gmres = false;
137 
138  // default GMRES parameters
139  m_mgmres = 100;
140  m_eps_gmres = 1.e-4;
141 
142  // some work space
143  m_spwork1.resize(m_nsp);
144  m_spwork2.resize(m_nsp);
145  m_spwork3.resize(m_nsp);
146 
147  // precompute and store log(epsilon_ij/k_B)
148  m_log_eps_k.resize(m_nsp, m_nsp);
149  // int j;
150  for (size_t i = 0; i < m_nsp; i++) {
151  for (size_t j = i; j < m_nsp; j++) {
152  m_log_eps_k(i,j) = log(tr.epsilon(i,j)/Boltzmann);
153  m_log_eps_k(j,i) = m_log_eps_k(i,j);
154  }
155  }
156 
157  // precompute and store constant parts of the Parker rotational
158  // collision number temperature correction
159  const doublereal sq298 = sqrt(298.0);
160  const doublereal kb298 = Boltzmann * 298.0;
161  m_sqrt_eps_k.resize(m_nsp);
162  for (size_t k = 0; k < m_nsp; k++) {
163  m_sqrt_eps_k[k] = sqrt(tr.eps[k]/Boltzmann);
164  m_frot_298[k] = Frot(tr.eps[k]/kb298,
165  m_sqrt_eps_k[k]/sq298);
166  }
167 
168  return true;
169 }
170 
171 //====================================================================================================================
172 
173 /****************** thermal conductivity **********************/
174 
175 /**
176  * @internal
177  */
179 {
181  doublereal sum = 0.0;
182  for (size_t k = 0; k < 2*m_nsp; k++) {
183  sum += m_b[k + m_nsp] * m_a[k + m_nsp];
184  }
185  return -4.0*sum;
186 }
187 //====================================================================================================================
188 // Return the thermal diffusion coefficients for the species
189 /*
190  *
191  * @param dt thermal diffusion coefficients
192  * (length = m_nsp)
193  */
194 void MultiTransport::getThermalDiffCoeffs(doublereal* const dt)
195 {
197  const doublereal c = 1.6/GasConstant;
198  for (size_t k = 0; k < m_nsp; k++) {
199  dt[k] = c * m_mw[k] * m_molefracs[k] * m_a[k];
200  }
201 }
202 //====================================================================================================================
203 
204 /**
205  * @internal
206  */
208 {
209  // if T has changed, update the temperature-dependent properties.
210  updateThermal_T();
211  update_C();
212 
213  // Copy the mole fractions twice into the last two blocks of
214  // the right-hand-side vector m_b. The first block of m_b was
215  // set to zero when it was created, and is not modified so
216  // doesn't need to be reset to zero.
217  for (size_t k = 0; k < m_nsp; k++) {
218  m_b[k] = 0.0;
219  m_b[k + m_nsp] = m_molefracs[k];
220  m_b[k + 2*m_nsp] = m_molefracs[k];
221  }
222 
223  // Set the right-hand side vector to zero in the 3rd block for
224  // all species with no internal energy modes. The
225  // corresponding third-block rows and columns will be set to
226  // zero, except on the diagonal of L01,01, where they are set
227  // to 1.0. This has the effect of eliminating these equations
228  // from the system, since the equation becomes: m_a[2*m_nsp +
229  // k] = 0.0.
230 
231  // Note that this differs from the Chemkin procedure, where
232  // all *monatomic* species are excluded. Since monatomic
233  // radicals can have non-zero internal heat capacities due to
234  // electronic excitation, they should be retained.
235  //
236  // But if CHEMKIN_COMPATIBILITY_MODE is defined, then all
237  // monatomic species are excluded.
238 
239  for (size_t k = 0; k < m_nsp; k++) {
240  if (!hasInternalModes(k)) {
241  m_b[2*m_nsp + k] = 0.0;
242  }
243  }
244 
245  // evaluate the submatrices of the L matrix
246  m_Lmatrix.resize(3*m_nsp, 3*m_nsp, 0.0);
247 
250  eval_L0001();
251  eval_L1000();
252  eval_L1010(DATA_PTR(m_molefracs));
253  eval_L1001(DATA_PTR(m_molefracs));
254  eval_L0100();
255  eval_L0110();
256  eval_L0101(DATA_PTR(m_molefracs));
257 
258  // Solve it using GMRES or LU decomposition. The last solution
259  // in m_a should provide a good starting guess, so convergence
260  // should be fast.
261 
262  //if (m_gmres) {
263  // gmres(m_mgmres, 3*m_nsp, m_Lmatrix, m_b.begin(),
264  // m_a.begin(), m_eps_gmres);
265  // m_lmatrix_soln_ok = true;
266  // m_l0000_ok = true; // L matrix not modified by GMRES
267  //}
268  //else {
269  copy(m_b.begin(), m_b.end(), m_a.begin());
270  try {
271  solve(m_Lmatrix, DATA_PTR(m_a));
272  } catch (CanteraError& err) {
273  err.save();
274  //if (info != 0) {
275  throw CanteraError("MultiTransport::solveLMatrixEquation",
276  "error in solving L matrix.");
277  }
278  m_lmatrix_soln_ok = true;
279  m_l0000_ok = false;
280  // L matrix is overwritten with LU decomposition
281  //}
282  m_lmatrix_soln_ok = true;
283 }
284 
285 //====================================================================================================================
286 // Get the species diffusive mass fluxes wrt to the mass averaged velocity,
287 // given the gradients in mole fraction and temperature
288 /*
289  * Units for the returned fluxes are kg m-2 s-1.
290  *
291  * @param ndim Number of dimensions in the flux expressions
292  * @param grad_T Gradient of the temperature
293  * (length = ndim)
294  * @param ldx Leading dimension of the grad_X array
295  * (usually equal to m_nsp but not always)
296  * @param grad_X Gradients of the mole fraction
297  * Flat vector with the m_nsp in the inner loop.
298  * length = ldx * ndim
299  * @param ldf Leading dimension of the fluxes array
300  * (usually equal to m_nsp but not always)
301  * @param fluxes Output of the diffusive mass fluxes
302  * Flat vector with the m_nsp in the inner loop.
303  * length = ldx * ndim
304  */
305 void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_T,
306  size_t ldx, const doublereal* const grad_X,
307  size_t ldf, doublereal* const fluxes)
308 {
309  // update the binary diffusion coefficients if necessary
310  update_T();
311  updateDiff_T();
312 
313  // If any component of grad_T is non-zero, then get the
314  // thermal diffusion coefficients
315  bool addThermalDiffusion = false;
316  for (size_t i = 0; i < ndim; i++) {
317  if (grad_T[i] != 0.0) {
318  addThermalDiffusion = true;
319  }
320  }
321  if (addThermalDiffusion) {
323  }
324 
325  const doublereal* y = m_thermo->massFractions();
326  doublereal rho = m_thermo->density();
327 
328  for (size_t i = 0; i < m_nsp; i++) {
329  double sum = 0.0;
330  for (size_t j = 0; j < m_nsp; j++) {
331  m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j);
332  sum += m_aa(i,j);
333  }
334  m_aa(i,i) -= sum;
335  }
336 
337  // enforce the condition \sum Y_k V_k = 0. This is done by replacing
338  // the flux equation with the largest gradx component in the first
339  // coordinate direction with the flux balance condition.
340  size_t jmax = 0;
341  doublereal gradmax = -1.0;
342  for (size_t j = 0; j < m_nsp; j++) {
343  if (fabs(grad_X[j]) > gradmax) {
344  gradmax = fabs(grad_X[j]);
345  jmax = j;
346  }
347  }
348 
349  // set the matrix elements in this row to the mass fractions,
350  // and set the entry in gradx to zero
351  for (size_t j = 0; j < m_nsp; j++) {
352  m_aa(jmax,j) = y[j];
353  }
354  vector_fp gsave(ndim), grx(ldx*m_nsp);
355  for (size_t n = 0; n < ldx*ndim; n++) {
356  grx[n] = grad_X[n];
357  }
358  //for (n = 0; n < ndim; n++) {
359  // gsave[n] = grad_X[jmax + n*ldx]; // save the input mole frac gradient
360  //grad_X[jmax + n*ldx] = 0.0;
361  // grx[jmax + n*ldx] = 0.0;
362  // }
363 
364  // copy grad_X to fluxes
365  const doublereal* gx;
366  for (size_t n = 0; n < ndim; n++) {
367  gx = grad_X + ldx*n;
368  copy(gx, gx + m_nsp, fluxes + ldf*n);
369  fluxes[jmax + n*ldf] = 0.0;
370  }
371 
372  // use LAPACK to solve the equations
373  int info=0;
374  ct_dgetrf(static_cast<int>(m_aa.nRows()),
375  static_cast<int>(m_aa.nColumns()), m_aa.ptrColumn(0),
376  static_cast<int>(m_aa.nRows()),
377  &m_aa.ipiv()[0], info);
378  if (info == 0) {
379  ct_dgetrs(ctlapack::NoTranspose,
380  static_cast<int>(m_aa.nRows()), ndim,
381  m_aa.ptrColumn(0), static_cast<int>(m_aa.nRows()),
382  &m_aa.ipiv()[0], fluxes, ldf, info);
383  if (info != 0) {
384  info += 100;
385  }
386  } else
387  throw CanteraError("MultiTransport::getSpeciesFluxes",
388  "Error in DGETRF");
389  if (info > 50)
390  throw CanteraError("MultiTransport::getSpeciesFluxes",
391  "Error in DGETRS");
392 
393  size_t offset;
394  doublereal pp = pressure_ig();
395 
396  // multiply diffusion velocities by rho * V to create
397  // mass fluxes, and restore the gradx elements that were
398  // modified
399  for (size_t n = 0; n < ndim; n++) {
400  offset = n*ldf;
401  for (size_t i = 0; i < m_nsp; i++) {
402  fluxes[i + offset] *= rho * y[i] / pp;
403  }
404  //grad_X[jmax + n*ldx] = gsave[n];
405  }
406 
407  // thermal diffusion
408  if (addThermalDiffusion) {
409  for (size_t n = 0; n < ndim; n++) {
410  offset = n*ldf;
411  doublereal grad_logt = grad_T[n]/m_temp;
412  for (size_t i = 0; i < m_nsp; i++) {
413  fluxes[i + offset] -= m_spwork[i]*grad_logt;
414  }
415  }
416  }
417 }
418 //====================================================================================================================
419 // Get the mass diffusional fluxes [kg/m^2/s] of the species, given the thermodynamic
420 // state at two nearby points.
421 /*
422  * The specific diffusional fluxes are calculated with reference to the mass averaged
423  * velocity. This is a one-dimensional vector
424  *
425  * @param state1 Array of temperature, density, and mass
426  * fractions for state 1.
427  * @param state2 Array of temperature, density, and mass
428  * fractions for state 2.
429  * @param delta Distance from state 1 to state 2 (m).
430  * @param fluxes Output mass fluxes of the species.
431  * (length = m_nsp)
432  */
433 void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* state2, doublereal delta,
434  doublereal* fluxes)
435 {
436  double* x1 = DATA_PTR(m_spwork1);
437  double* x2 = DATA_PTR(m_spwork2);
438  double* x3 = DATA_PTR(m_spwork3);
439  size_t n, nsp = m_thermo->nSpecies();
440  m_thermo->restoreState(nsp+2, state1);
441  double p1 = m_thermo->pressure();
442  double t1 = state1[0];
444 
445  m_thermo->restoreState(nsp+2, state2);
446  double p2 = m_thermo->pressure();
447  double t2 = state2[0];
449 
450  double p = 0.5*(p1 + p2);
451  double t = 0.5*(state1[0] + state2[0]);
452 
453  for (n = 0; n < nsp; n++) {
454  x3[n] = 0.5*(x1[n] + x2[n]);
455  }
456  m_thermo->setState_TPX(t, p, x3);
458 
459  // update the binary diffusion coefficients if necessary
460  update_T();
461  updateDiff_T();
462 
463  // If there is a temperature gradient, then get the
464  // thermal diffusion coefficients
465 
466  bool addThermalDiffusion = false;
467  if (state1[0] != state2[0]) {
468  addThermalDiffusion = true;
470  }
471 
472  const doublereal* y = m_thermo->massFractions();
473  doublereal rho = m_thermo->density();
474 
475  for (size_t i = 0; i < m_nsp; i++) {
476  doublereal sum = 0.0;
477  for (size_t j = 0; j < m_nsp; j++) {
478  m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j);
479  sum += m_aa(i,j);
480  }
481  m_aa(i,i) -= sum;
482  }
483 
484  // enforce the condition \sum Y_k V_k = 0. This is done by
485  // replacing the flux equation with the largest gradx
486  // component with the flux balance condition.
487  size_t jmax = 0;
488  doublereal gradmax = -1.0;
489  for (size_t j = 0; j < m_nsp; j++) {
490  if (fabs(x2[j] - x1[j]) > gradmax) {
491  gradmax = fabs(x1[j] - x2[j]);
492  jmax = j;
493  }
494  }
495 
496  // set the matrix elements in this row to the mass fractions,
497  // and set the entry in gradx to zero
498 
499  for (size_t j = 0; j < m_nsp; j++) {
500  m_aa(jmax,j) = y[j];
501  fluxes[j] = x2[j] - x1[j];
502  }
503  fluxes[jmax] = 0.0;
504 
505  // use LAPACK to solve the equations
506  int info=0;
507  size_t nr = m_aa.nRows();
508  size_t nc = m_aa.nColumns();
509 
510  ct_dgetrf(nr, nc, m_aa.ptrColumn(0), nr, &m_aa.ipiv()[0], info);
511  if (info == 0) {
512  int ndim = 1;
513  ct_dgetrs(ctlapack::NoTranspose, nr, ndim,
514  m_aa.ptrColumn(0), nr, &m_aa.ipiv()[0], fluxes, nr, info);
515  if (info != 0)
516  throw CanteraError("MultiTransport::getMassFluxes",
517  "Error in DGETRS. Info = "+int2str(info));
518  } else
519  throw CanteraError("MultiTransport::getMassFluxes",
520  "Error in DGETRF. Info = "+int2str(info));
521 
522  doublereal pp = pressure_ig();
523 
524  // multiply diffusion velocities by rho * Y_k to create
525  // mass fluxes, and divide by pressure
526  for (size_t i = 0; i < m_nsp; i++) {
527  fluxes[i] *= rho * y[i] / pp;
528  }
529 
530  // thermal diffusion
531  if (addThermalDiffusion) {
532  doublereal grad_logt = (t2 - t1)/m_temp;
533  for (size_t i = 0; i < m_nsp; i++) {
534  fluxes[i] -= m_spwork[i]*grad_logt;
535  }
536  }
537 }
538 //====================================================================================================================
539 void MultiTransport::getMolarFluxes(const doublereal* const state1,
540  const doublereal* const state2,
541  const doublereal delta,
542  doublereal* const fluxes)
543 {
544  getMassFluxes(state1, state2, delta, fluxes);
545  for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
546  fluxes[k] /= m_mw[k];
547  }
548 }
549 //====================================================================================================================
550 // Set the solution method for inverting the L matrix
551 /*
552  * @param method enum TRANSOLVE_TYPE Either use direct or TRANSOLVE_GMRES
553  */
555 {
556  if (method == TRANSOLVE_GMRES) {
557  m_gmres = true;
558  } else {
559  m_gmres = false;
560  }
561 }
562 //====================================================================================================================
563 void MultiTransport::setOptions_GMRES(int m, doublereal eps)
564 {
565  if (m > 0) {
566  m_mgmres = m;
567  }
568  if (eps > 0.0) {
569  m_eps_gmres = eps;
570  }
571 }
572 //====================================================================================================================
573 void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
574 {
575  doublereal p = pressure_ig();
576 
577  // update the mole fractions
578  update_C();
579 
580  // update the binary diffusion coefficients
581  update_T();
582  updateThermal_T();
583 
584  // evaluate L0000 if the temperature or concentrations have
585  // changed since it was last evaluated.
586  if (!m_l0000_ok) {
588  }
589 
590  // invert L00,00
591  int ierr = invert(m_Lmatrix, m_nsp);
592  if (ierr != 0) {
593  throw CanteraError("MultiTransport::getMultiDiffCoeffs",
594  string(" invert returned ierr = ")+int2str(ierr));
595  }
596  m_l0000_ok = false; // matrix is overwritten by inverse
597 
598  //doublereal pres = m_thermo->pressure();
599  doublereal prefactor = 16.0 * m_temp
600  * m_thermo->meanMolecularWeight()/(25.0 * p);
601  doublereal c;
602 
603  for (size_t i = 0; i < m_nsp; i++) {
604  for (size_t j = 0; j < m_nsp; j++) {
605  c = prefactor/m_mw[j];
606  d[ld*j + i] = c*m_molefracs[i]*
607  (m_Lmatrix(i,j) - m_Lmatrix(i,i));
608  }
609  }
610 }
611 //====================================================================================================================
612 
613 
615 {
616  if (m_temp == m_thermo->temperature()) {
617  return;
618  }
619 
620  GasTransport::update_T();
621 
622  // temperature has changed, so polynomial fits will need to be
623  // redone, and the L matrix reevaluated.
624  m_abc_ok = false;
625  m_lmatrix_soln_ok = false;
626  m_l0000_ok = false;
627 }
628 
630 {
631  // signal that concentration-dependent quantities will need to
632  // be recomputed before use, and update the local mole
633  // fraction array.
634  m_l0000_ok = false;
635  m_lmatrix_soln_ok = false;
637 
638  // add an offset to avoid a pure species condition
639  // (check - this may be unnecessary)
640  for (size_t k = 0; k < m_nsp; k++) {
642  }
643 }
644 
645 /*************************************************************************
646  *
647  * methods to update temperature-dependent properties
648  *
649  *************************************************************************/
650 
652 {
653  if (m_thermal_tlast == m_thermo->temperature()) {
654  return;
655  }
656  // we need species viscosities and binary diffusion coefficients
658  updateDiff_T();
659 
660  // evaluate polynomial fits for A*, B*, C*
661  doublereal z;
662  int ipoly;
663  for (size_t i = 0; i < m_nsp; i++) {
664  for (size_t j = i; j < m_nsp; j++) {
665  z = m_logt - m_log_eps_k(i,j);
666  ipoly = m_poly[i][j];
667  if (m_mode == CK_Mode) {
668  m_om22(i,j) = poly6(z, DATA_PTR(m_om22_poly[ipoly]));
669  m_astar(i,j) = poly6(z, DATA_PTR(m_astar_poly[ipoly]));
670  m_bstar(i,j) = poly6(z, DATA_PTR(m_bstar_poly[ipoly]));
671  m_cstar(i,j) = poly6(z, DATA_PTR(m_cstar_poly[ipoly]));
672  } else {
673  m_om22(i,j) = poly8(z, DATA_PTR(m_om22_poly[ipoly]));
674  m_astar(i,j) = poly8(z, DATA_PTR(m_astar_poly[ipoly]));
675  m_bstar(i,j) = poly8(z, DATA_PTR(m_bstar_poly[ipoly]));
676  m_cstar(i,j) = poly8(z, DATA_PTR(m_cstar_poly[ipoly]));
677  }
678  m_om22(j,i) = m_om22(i,j);
679  m_astar(j,i) = m_astar(i,j);
680  m_bstar(j,i) = m_bstar(i,j);
681  m_cstar(j,i) = m_cstar(i,j);
682  }
683  }
684  m_abc_ok = true;
685 
686  // evaluate the temperature-dependent rotational relaxation rate
687  doublereal tr, sqtr;
688  for (size_t k = 0; k < m_nsp; k++) {
689  tr = m_eps[k]/ m_kbt;
690  sqtr = m_sqrt_eps_k[k] / m_sqrt_t;
691  m_rotrelax[k] = std::max(1.0,m_zrot[k]) * m_frot_298[k]/Frot(tr, sqtr);
692  }
693 
694  doublereal d;
695  doublereal c = 1.2*GasConstant*m_temp;
696  for (size_t k = 0; k < m_nsp; k++) {
697  d = c * m_visc[k] * m_astar(k,k)/m_mw[k];
698  m_bdiff(k,k) = d;
699  }
700 
701  // Calculate the internal heat capacities by subtracting off the translational contributions
702  /*
703  * HKM Exploratory comment:
704  * The translational component is 1.5
705  * The rotational component is 1.0 for a linear molecule and 1.5 for a nonlinear molecule
706  * and zero for a monatomic.
707  * Chemkin has traditionally subtracted 1.5 here (SAND86-8246).
708  * The original Dixon-Lewis paper subtracted 1.5 here.
709  */
710  const vector_fp& cp = ((IdealGasPhase*)m_thermo)->cp_R_ref();
711  for (size_t k = 0; k < m_nsp; k++) {
712  m_cinternal[k] = cp[k] - 2.5;
713  }
714 
715  // m_thermo->update_T(m_update_thermal_T);
716  m_thermal_tlast = m_thermo->temperature();
717 }
718 
719 //====================================================================================================================
720 /*
721  * This function returns a Transport data object for a given species.
722  *
723  */
725 getGasTransportData(int kSpecies) {
726  struct GasTransportData td;
727  td.speciesName = m_thermo->speciesName(kSpecies);
728 
729  td.geometry = 2;
730  if (m_crot[kSpecies] == 0.0) {
731  td.geometry = 0;
732  } else if (m_crot[kSpecies] == 1.0) {
733  td.geometry = 1;
734  }
735  td.wellDepth = m_eps[kSpecies] / Boltzmann;
736  td.dipoleMoment = m_dipoleDiag[kSpecies] * 1.0E25 / SqrtTen;
737  td.diameter = m_diam(kSpecies, kSpecies) * 1.0E10;
738  td.polarizability = m_alpha[kSpecies] * 1.0E30;
739  td.rotRelaxNumber = m_zrot[kSpecies];
740 
741  return td;
742 }
743 //====================================================================================================================
744 }