Cantera  2.2.1
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
13
14 using namespace std;
15
16 namespace Cantera
17 {
18
19 ///////////////////// helper functions /////////////////////////
20
21 /**
22  * The Parker temperature correction to the rotational collision number.
23  *
24  * @param tr Reduced temperature \f$\epsilon/kT \f$
25  * @param sqtr square root of tr.
26  */
27 inline doublereal Frot(doublereal tr, doublereal sqtr)
28 {
29  const doublereal c1 = 0.5*sqrt(Pi)*Pi;
30  const doublereal c2 = 0.25*Pi*Pi + 2.0;
31  const doublereal c3 = sqrt(Pi)*Pi;
32  return 1.0 + c1*sqtr + c2*tr + c3*sqtr*tr;
33 }
34
35 //////////////////// class MultiTransport methods //////////////
36
37 MultiTransport::MultiTransport(thermo_t* thermo)
38  : GasTransport(thermo)
39 {
40 }
41
42 void MultiTransport::init(ThermoPhase* thermo, int mode, int log_level)
43 {
44  GasTransport::init(thermo, mode, log_level);
45
46  // the L matrix
47  m_Lmatrix.resize(3*m_nsp, 3*m_nsp);
48  m_a.resize(3*m_nsp, 1.0);
49  m_b.resize(3*m_nsp, 0.0);
50  m_aa.resize(m_nsp, m_nsp, 0.0);
51  m_molefracs_last.resize(m_nsp, -1.0);
52
53  m_frot_298.resize(m_nsp);
54  m_rotrelax.resize(m_nsp);
55
56  m_cinternal.resize(m_nsp);
57
62
63  // set flags all false
64  m_abc_ok = false;
65  m_l0000_ok = false;
66  m_lmatrix_soln_ok = false;
67
68  m_thermal_tlast = 0.0;
69
70  // some work space
71  m_spwork1.resize(m_nsp);
72  m_spwork2.resize(m_nsp);
73  m_spwork3.resize(m_nsp);
74
75  // precompute and store log(epsilon_ij/k_B)
76  m_log_eps_k.resize(m_nsp, m_nsp);
77  // int j;
78  for (size_t i = 0; i < m_nsp; i++) {
79  for (size_t j = i; j < m_nsp; j++) {
80  m_log_eps_k(i,j) = log(m_epsilon(i,j)/Boltzmann);
81  m_log_eps_k(j,i) = m_log_eps_k(i,j);
82  }
83  }
84
85  // precompute and store constant parts of the Parker rotational
86  // collision number temperature correction
87  const doublereal sq298 = sqrt(298.0);
88  const doublereal kb298 = Boltzmann * 298.0;
89  m_sqrt_eps_k.resize(m_nsp);
90  for (size_t k = 0; k < m_nsp; k++) {
91  m_sqrt_eps_k[k] = sqrt(m_eps[k]/Boltzmann);
92  m_frot_298[k] = Frot(m_eps[k]/kb298, m_sqrt_eps_k[k]/sq298);
93  }
94 }
95
97 {
98  solveLMatrixEquation();
99  doublereal sum = 0.0;
100  for (size_t k = 0; k < 2*m_nsp; k++) {
101  sum += m_b[k + m_nsp] * m_a[k + m_nsp];
102  }
103  return -4.0*sum;
104 }
105
106 void MultiTransport::getThermalDiffCoeffs(doublereal* const dt)
107 {
108  solveLMatrixEquation();
109  const doublereal c = 1.6/GasConstant;
110  for (size_t k = 0; k < m_nsp; k++) {
111  dt[k] = c * m_mw[k] * m_molefracs[k] * m_a[k];
112  }
113 }
114
115 void MultiTransport::solveLMatrixEquation()
116 {
117  // if T has changed, update the temperature-dependent properties.
118  updateThermal_T();
119  update_C();
120  if (m_lmatrix_soln_ok) {
121  return;
122  }
123
124  // Copy the mole fractions twice into the last two blocks of
125  // the right-hand-side vector m_b. The first block of m_b was
126  // set to zero when it was created, and is not modified so
127  // doesn't need to be reset to zero.
128  for (size_t k = 0; k < m_nsp; k++) {
129  m_b[k] = 0.0;
130  m_b[k + m_nsp] = m_molefracs[k];
131  m_b[k + 2*m_nsp] = m_molefracs[k];
132  }
133
134  // Set the right-hand side vector to zero in the 3rd block for
135  // all species with no internal energy modes. The
136  // corresponding third-block rows and columns will be set to
137  // zero, except on the diagonal of L01,01, where they are set
138  // to 1.0. This has the effect of eliminating these equations
139  // from the system, since the equation becomes: m_a[2*m_nsp +
140  // k] = 0.0.
141
142  // Note that this differs from the Chemkin procedure, where
143  // all *monatomic* species are excluded. Since monatomic
144  // radicals can have non-zero internal heat capacities due to
145  // electronic excitation, they should be retained.
146
147  for (size_t k = 0; k < m_nsp; k++) {
148  if (!hasInternalModes(k)) {
149  m_b[2*m_nsp + k] = 0.0;
150  }
151  }
152
153  // evaluate the submatrices of the L matrix
154  m_Lmatrix.resize(3*m_nsp, 3*m_nsp, 0.0);
155
156  //! Evaluate the upper-left block of the L matrix.
159  eval_L0001();
160  eval_L1000();
161  eval_L1010(DATA_PTR(m_molefracs));
162  eval_L1001(DATA_PTR(m_molefracs));
163  eval_L0100();
164  eval_L0110();
165  eval_L0101(DATA_PTR(m_molefracs));
166
167  // Solve it using GMRES or LU decomposition. The last solution
168  // in m_a should provide a good starting guess, so convergence
169  // should be fast.
170
171  copy(m_b.begin(), m_b.end(), m_a.begin());
172  try {
173  solve(m_Lmatrix, DATA_PTR(m_a));
174  } catch (CanteraError& err) {
175  err.save();
176  throw CanteraError("MultiTransport::solveLMatrixEquation",
177  "error in solving L matrix.");
178  }
179  m_lmatrix_soln_ok = true;
181  // L matrix is overwritten with LU decomposition
182  m_l0000_ok = false;
183 }
184
185 void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_T,
186  size_t ldx, const doublereal* const grad_X,
187  size_t ldf, doublereal* const fluxes)
188 {
189  // update the binary diffusion coefficients if necessary
190  update_T();
191  updateDiff_T();
192
193  // If any component of grad_T is non-zero, then get the
194  // thermal diffusion coefficients
196  for (size_t i = 0; i < ndim; i++) {
197  if (grad_T[i] != 0.0) {
199  }
200  }
203  }
204
205  const doublereal* y = m_thermo->massFractions();
206  doublereal rho = m_thermo->density();
207
208  for (size_t i = 0; i < m_nsp; i++) {
209  double sum = 0.0;
210  for (size_t j = 0; j < m_nsp; j++) {
211  m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j);
212  sum += m_aa(i,j);
213  }
214  m_aa(i,i) -= sum;
215  }
216
217  // enforce the condition \sum Y_k V_k = 0. This is done by replacing
218  // the flux equation with the largest gradx component in the first
219  // coordinate direction with the flux balance condition.
220  size_t jmax = 0;
222  for (size_t j = 0; j < m_nsp; j++) {
225  jmax = j;
226  }
227  }
228
229  // set the matrix elements in this row to the mass fractions,
230  // and set the entry in gradx to zero
231  for (size_t j = 0; j < m_nsp; j++) {
232  m_aa(jmax,j) = y[j];
233  }
234  vector_fp gsave(ndim), grx(ldx*m_nsp);
235  for (size_t n = 0; n < ldx*ndim; n++) {
237  }
238
239  // copy grad_X to fluxes
240  const doublereal* gx;
241  for (size_t n = 0; n < ndim; n++) {
242  gx = grad_X + ldx*n;
243  copy(gx, gx + m_nsp, fluxes + ldf*n);
244  fluxes[jmax + n*ldf] = 0.0;
245  }
246
247  // use LAPACK to solve the equations
248  int info = m_aa.factor();
249  if (info) {
250  throw CanteraError("MultiTransport::getSpeciesFluxes",
251  "Error factorizing matrix.");
252  }
253  info = m_aa.solve(fluxes, ndim, ldf);
254  if (info) {
255  throw CanteraError("MultiTransport::getSpeciesFluxes",
256  "Error solving linear system.");
257  }
258
259  size_t offset;
260  doublereal pp = pressure_ig();
261
262  // multiply diffusion velocities by rho * V to create
263  // mass fluxes, and restore the gradx elements that were
264  // modified
265  for (size_t n = 0; n < ndim; n++) {
266  offset = n*ldf;
267  for (size_t i = 0; i < m_nsp; i++) {
268  fluxes[i + offset] *= rho * y[i] / pp;
269  }
270  }
271
272  // thermal diffusion
274  for (size_t n = 0; n < ndim; n++) {
275  offset = n*ldf;
277  for (size_t i = 0; i < m_nsp; i++) {
278  fluxes[i + offset] -= m_spwork[i]*grad_logt;
279  }
280  }
281  }
282 }
283
284 void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* state2, doublereal delta,
285  doublereal* fluxes)
286 {
287  double* x1 = DATA_PTR(m_spwork1);
288  double* x2 = DATA_PTR(m_spwork2);
289  double* x3 = DATA_PTR(m_spwork3);
290  size_t n, nsp = m_thermo->nSpecies();
291  m_thermo->restoreState(nsp+2, state1);
292  double p1 = m_thermo->pressure();
293  double t1 = state1[0];
295
296  m_thermo->restoreState(nsp+2, state2);
297  double p2 = m_thermo->pressure();
298  double t2 = state2[0];
300
301  double p = 0.5*(p1 + p2);
302  double t = 0.5*(state1[0] + state2[0]);
303
304  for (n = 0; n < nsp; n++) {
305  x3[n] = 0.5*(x1[n] + x2[n]);
306  }
307  m_thermo->setState_TPX(t, p, x3);
309
310  // update the binary diffusion coefficients if necessary
311  update_T();
312  updateDiff_T();
313
314  // If there is a temperature gradient, then get the
315  // thermal diffusion coefficients
316
318  if (state1[0] != state2[0]) {
321  }
322
323  const doublereal* y = m_thermo->massFractions();
324  doublereal rho = m_thermo->density();
325
326  for (size_t i = 0; i < m_nsp; i++) {
327  doublereal sum = 0.0;
328  for (size_t j = 0; j < m_nsp; j++) {
329  m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j);
330  sum += m_aa(i,j);
331  }
332  m_aa(i,i) -= sum;
333  }
334
335  // enforce the condition \sum Y_k V_k = 0. This is done by
336  // replacing the flux equation with the largest gradx
337  // component with the flux balance condition.
338  size_t jmax = 0;
340  for (size_t j = 0; j < m_nsp; j++) {
341  if (fabs(x2[j] - x1[j]) > gradmax) {
342  gradmax = fabs(x1[j] - x2[j]);
343  jmax = j;
344  }
345  }
346
347  // set the matrix elements in this row to the mass fractions,
348  // and set the entry in gradx to zero
349
350  for (size_t j = 0; j < m_nsp; j++) {
351  m_aa(jmax,j) = y[j];
352  fluxes[j] = x2[j] - x1[j];
353  }
354  fluxes[jmax] = 0.0;
355
356  // Solve the equations
357  int info = m_aa.factor();
358  if (info) {
359  throw CanteraError("MultiTransport::getMassFluxes",
360  "Error in factorization. Info = "+int2str(info));
361  }
362  info = m_aa.solve(fluxes);
363  if (info) {
364  throw CanteraError("MultiTransport::getMassFluxes",
365  "Error in linear solve. Info = "+int2str(info));
366  }
367
368  doublereal pp = pressure_ig();
369
370  // multiply diffusion velocities by rho * Y_k to create
371  // mass fluxes, and divide by pressure
372  for (size_t i = 0; i < m_nsp; i++) {
373  fluxes[i] *= rho * y[i] / pp;
374  }
375
376  // thermal diffusion
378  doublereal grad_logt = (t2 - t1)/m_temp;
379  for (size_t i = 0; i < m_nsp; i++) {
381  }
382  }
383 }
384
385 void MultiTransport::getMolarFluxes(const doublereal* const state1,
386  const doublereal* const state2,
387  const doublereal delta,
388  doublereal* const fluxes)
389 {
390  getMassFluxes(state1, state2, delta, fluxes);
391  for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
392  fluxes[k] /= m_mw[k];
393  }
394 }
395
396 void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
397 {
398  doublereal p = pressure_ig();
399
400  // update the mole fractions
401  update_C();
402
403  // update the binary diffusion coefficients
404  update_T();
405  updateThermal_T();
406
407  // evaluate L0000 if the temperature or concentrations have
408  // changed since it was last evaluated.
409  if (!m_l0000_ok) {
411  }
412
413  // invert L00,00
414  int ierr = invert(m_Lmatrix, m_nsp);
415  if (ierr != 0) {
416  throw CanteraError("MultiTransport::getMultiDiffCoeffs",
417  string(" invert returned ierr = ")+int2str(ierr));
418  }
419  m_l0000_ok = false; // matrix is overwritten by inverse
420  m_lmatrix_soln_ok = false;
421
422  doublereal prefactor = 16.0 * m_temp
423  * m_thermo->meanMolecularWeight()/(25.0 * p);
424  doublereal c;
425
426  for (size_t i = 0; i < m_nsp; i++) {
427  for (size_t j = 0; j < m_nsp; j++) {
428  c = prefactor/m_mw[j];
429  d[ld*j + i] = c*m_molefracs[i]*
430  (m_Lmatrix(i,j) - m_Lmatrix(i,i));
431  }
432  }
433 }
434
436 {
437  if (m_temp == m_thermo->temperature()) {
438  return;
439  }
440
441  GasTransport::update_T();
442
443  // temperature has changed, so polynomial fits will need to be
444  // redone, and the L matrix reevaluated.
445  m_abc_ok = false;
446  m_lmatrix_soln_ok = false;
447  m_l0000_ok = false;
448 }
449
451 {
452  // Update the local mole fraction array
454
455  for (size_t k = 0; k < m_nsp; k++) {
456  // add an offset to avoid a pure species condition
457  m_molefracs[k] = std::max(Tiny, m_molefracs[k]);
458  if (m_molefracs[k] != m_molefracs_last[k]) {
459  // If any mole fractions have changed, signal that concentration-
460  // dependent quantities will need to be recomputed before use.
461  m_l0000_ok = false;
462  m_lmatrix_soln_ok = false;
463  }
464  }
465 }
466
468 {
469  if (m_thermal_tlast == m_thermo->temperature()) {
470  return;
471  }
472  // we need species viscosities and binary diffusion coefficients
474  updateDiff_T();
475
476  // evaluate polynomial fits for A*, B*, C*
477  doublereal z;
478  int ipoly;
479  for (size_t i = 0; i < m_nsp; i++) {
480  for (size_t j = i; j < m_nsp; j++) {
481  z = m_logt - m_log_eps_k(i,j);
482  ipoly = m_poly[i][j];
483  if (m_mode == CK_Mode) {
484  m_om22(i,j) = poly6(z, DATA_PTR(m_omega22_poly[ipoly]));
485  m_astar(i,j) = poly6(z, DATA_PTR(m_astar_poly[ipoly]));
486  m_bstar(i,j) = poly6(z, DATA_PTR(m_bstar_poly[ipoly]));
487  m_cstar(i,j) = poly6(z, DATA_PTR(m_cstar_poly[ipoly]));
488  } else {
489  m_om22(i,j) = poly8(z, DATA_PTR(m_omega22_poly[ipoly]));
490  m_astar(i,j) = poly8(z, DATA_PTR(m_astar_poly[ipoly]));
491  m_bstar(i,j) = poly8(z, DATA_PTR(m_bstar_poly[ipoly]));
492  m_cstar(i,j) = poly8(z, DATA_PTR(m_cstar_poly[ipoly]));
493  }
494  m_om22(j,i) = m_om22(i,j);
495  m_astar(j,i) = m_astar(i,j);
496  m_bstar(j,i) = m_bstar(i,j);
497  m_cstar(j,i) = m_cstar(i,j);
498  }
499  }
500  m_abc_ok = true;
501
502  // evaluate the temperature-dependent rotational relaxation rate
503  doublereal tr, sqtr;
504  for (size_t k = 0; k < m_nsp; k++) {
505  tr = m_eps[k]/ m_kbt;
506  sqtr = m_sqrt_eps_k[k] / m_sqrt_t;
507  m_rotrelax[k] = std::max(1.0,m_zrot[k]) * m_frot_298[k]/Frot(tr, sqtr);
508  }
509
510  doublereal d;
511  doublereal c = 1.2*GasConstant*m_temp;
512  for (size_t k = 0; k < m_nsp; k++) {
513  d = c * m_visc[k] * m_astar(k,k)/m_mw[k];
514  m_bdiff(k,k) = d;
515  }
516
517  // Calculate the internal heat capacities by subtracting off the translational contributions
518  /*
519  * HKM Exploratory comment:
520  * The translational component is 1.5
521  * The rotational component is 1.0 for a linear molecule and 1.5 for a nonlinear molecule
522  * and zero for a monatomic.
523  * Chemkin has traditionally subtracted 1.5 here (SAND86-8246).
524  * The original Dixon-Lewis paper subtracted 1.5 here.
525  */
526  vector_fp cp(m_thermo->nSpecies());
527  m_thermo->getCp_R_ref(&cp[0]);
528
529  for (size_t k = 0; k < m_nsp; k++) {
530  m_cinternal[k] = cp[k] - 2.5;
531  }
532
533  m_thermal_tlast = m_thermo->temperature();
534 }
535
536 //! Constant to compare dimensionless heat capacities against zero
537 static const doublereal Min_C_Internal = 0.001;
538
539 bool MultiTransport::hasInternalModes(size_t j)
540 {
541  return (m_cinternal[j] > Min_C_Internal);
542 }
543
544 void MultiTransport::eval_L0000(const doublereal* const x)
545 {
546  doublereal prefactor = 16.0*m_temp/25.0;
547  doublereal sum;
548  for (size_t i = 0; i < m_nsp; i++) {
549  // subtract-off the k=i term to account for the first delta
550  // function in Eq. (12.121)
551
552  sum = -x[i]/m_bdiff(i,i);
553  for (size_t k = 0; k < m_nsp; k++) {
554  sum += x[k]/m_bdiff(i,k);
555  }
556
557  sum /= m_mw[i];
558  for (size_t j = 0; j != m_nsp; ++j) {
559  m_Lmatrix(i,j) = prefactor * x[j]
560  * (m_mw[j] * sum + x[i]/m_bdiff(i,j));
561  }
562  // diagonal term is zero
563  m_Lmatrix(i,i) = 0.0;
564  }
565 }
566
567 void MultiTransport::eval_L0010(const doublereal* const x)
568 {
569  doublereal prefactor = 1.6*m_temp;
570
571  doublereal sum, wj, xj;
572  for (size_t j = 0; j < m_nsp; j++) {
573  xj = x[j];
574  wj = m_mw[j];
575  sum = 0.0;
576  for (size_t i = 0; i < m_nsp; i++) {
577  m_Lmatrix(i,j + m_nsp) = - prefactor * x[i] * xj * m_mw[i] *
578  (1.2 * m_cstar(j,i) - 1.0) /
579  ((wj + m_mw[i]) * m_bdiff(j,i));
580
581  // the next term is independent of "j";
582  // need to do it for the "j,j" term
583  sum -= m_Lmatrix(i,j+m_nsp);
584  }
585  m_Lmatrix(j,j+m_nsp) += sum;
586  }
587 }
588
590 {
591  for (size_t j = 0; j < m_nsp; j++) {
592  for (size_t i = 0; i < m_nsp; i++) {
593  m_Lmatrix(i+m_nsp,j) = m_Lmatrix(j,i+m_nsp);
594  }
595  }
596 }
597
598 void MultiTransport::eval_L1010(const doublereal* x)
599 {
600  const doublereal fiveover3pi = 5.0/(3.0*Pi);
601  doublereal prefactor = (16.0*m_temp)/25.0;
602
603  doublereal constant1, wjsq, constant2, constant3, constant4,
604  fourmj, threemjsq, sum, sumwij;;
605  doublereal term1, term2;
606
607  for (size_t j = 0; j < m_nsp; j++) {
608
609  // get constant terms that depend on just species "j"
610
611  constant1 = prefactor*x[j];
612  wjsq = m_mw[j]*m_mw[j];
613  constant2 = 13.75*wjsq;
614  constant3 = m_crot[j]/m_rotrelax[j];
615  constant4 = 7.5*wjsq;
616  fourmj = 4.0*m_mw[j];
617  threemjsq = 3.0*m_mw[j]*m_mw[j];
618  sum = 0.0;
619  for (size_t i = 0; i < m_nsp; i++) {
620
621  sumwij = m_mw[i] + m_mw[j];
622  term1 = m_bdiff(i,j) * sumwij*sumwij;
623  term2 = fourmj*m_astar(i,j)*(1.0 + fiveover3pi*
624  (constant3 +
625  (m_crot[i]/m_rotrelax[i]))); // see Eq. (12.125)
626
627  m_Lmatrix(i+m_nsp,j+m_nsp) = constant1*x[i]*m_mw[i] /(m_mw[j]*term1) *
628  (constant2 - threemjsq*m_bstar(i,j)
629  - term2*m_mw[j]);
630
631  sum += x[i] /(term1) *
632  (constant4 + m_mw[i]*m_mw[i]*
633  (6.25 - 3.0*m_bstar(i,j)) + term2*m_mw[i]);
634  }
635
636  m_Lmatrix(j+m_nsp,j+m_nsp) -= sum*constant1;
637  }
638 }
639
640 void MultiTransport::eval_L1001(const doublereal* x)
641 {
642  doublereal prefactor = 32.00*m_temp/(5.00*Pi);
643  doublereal constant, sum;
644  size_t n2 = 2*m_nsp;
645  int npoly = 0;
646  for (size_t j = 0; j < m_nsp; j++) {
647  // collect terms that depend only on "j"
648  if (hasInternalModes(j)) {
649  constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]);
650  sum = 0.0;
651  for (size_t i = 0; i < m_nsp; i++) {
652  // see Eq. (12.127)
653  m_Lmatrix(i+m_nsp,j+n2) = constant * m_astar(j,i) * x[i] /
654  ((m_mw[j] + m_mw[i]) * m_bdiff(j,i));
655  sum += m_Lmatrix(i+m_nsp,j+n2);
656  }
657  npoly++;
658  m_Lmatrix(j+m_nsp,j+n2) += sum;
659  } else {
660  for (size_t i = 0; i < m_nsp; i++) {
661  m_Lmatrix(i+m_nsp,j+n2) = 0.0;
662  }
663  }
664  }
665 }
666
667 void MultiTransport::eval_L0001()
668 {
669  size_t n2 = 2*m_nsp;
670  for (size_t j = 0; j < m_nsp; j++) {
671  for (size_t i = 0; i < m_nsp; i++) {
672  m_Lmatrix(i,j+n2) = 0.0;
673  }
674  }
675 }
676
677 void MultiTransport::eval_L0100()
678 {
679  size_t n2 = 2*m_nsp;
680  for (size_t j = 0; j < m_nsp; j++)
681  for (size_t i = 0; i < m_nsp; i++) {
682  m_Lmatrix(i+n2,j) = 0.0; // see Eq. (12.123)
683  }
684 }
685
686 void MultiTransport::eval_L0110()
687 {
688  size_t n2 = 2*m_nsp;
689  for (size_t j = 0; j < m_nsp; j++)
690  for (size_t i = 0; i < m_nsp; i++) {
691  m_Lmatrix(i+n2,j+m_nsp) = m_Lmatrix(j+m_nsp,i+n2); // see Eq. (12.123)
692  }
693 }
694
695 void MultiTransport::eval_L0101(const doublereal* x)
696 {
697  const doublereal fivepi = 5.00*Pi;
698  const doublereal eightoverpi = 8.0 / Pi;
699
700  doublereal prefactor = 4.00*m_temp;
701  size_t n2 = 2*m_nsp;
702  doublereal constant1, constant2, diff_int, sum;
703  for (size_t i = 0; i < m_nsp; i++) {
704  if (hasInternalModes(i)) {
705  // collect terms that depend only on "i"
706  constant1 = prefactor*x[i]/m_cinternal[i];
707  constant2 = 12.00*m_mw[i]*m_crot[i] /
708  (fivepi*m_cinternal[i]*m_rotrelax[i]);
709  sum = 0.0;
710  for (size_t k = 0; k < m_nsp; k++) {
711  // see Eq. (12.131)
712  diff_int = m_bdiff(i,k);
713  m_Lmatrix(k+n2,i+n2) = 0.0;
714  sum += x[k]/diff_int;
715  if (k != i) sum += x[k]*m_astar(i,k)*constant2 /
716  (m_mw[k]*diff_int);
717  }
718  // see Eq. (12.130)
719  m_Lmatrix(i+n2,i+n2) =
720  - eightoverpi*m_mw[i]*x[i]*x[i]*m_crot[i] /
721  (m_cinternal[i]*m_cinternal[i]*GasConstant*m_visc[i]*m_rotrelax[i])
722  - constant1*sum;
723  } else {
724  for (size_t k = 0; k < m_nsp; k++) {
725  m_Lmatrix(i+n2,i+n2) = 1.0;
726  }
727  }
728  }
729 }
730
731 }
virtual doublereal thermalConductivity()
Returns the mixture thermal conductivity in W/m/K.
void eval_L1000()
Evaluate the L1000 matrices.
Update the pure-species viscosities.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:608
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
DenseMatrix m_bstar
Dense matrix for bstar.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
void restoreState(const vector_fp &state)
Restore a state saved on a previous call to saveState.
Definition: Phase.cpp:314
virtual void getThermalDiffCoeffs(doublereal *const dt)
Return the thermal diffusion coefficients (kg/m/s)
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
Interface for class MultiTransport.
DenseMatrix m_astar
Dense matrix for astar.
virtual void getMassFluxes(const doublereal *state1, const doublereal *state2, doublereal delta, doublereal *fluxes)
Get the mass diffusional fluxes [kg/m^2/s] of the species, given the thermodynamic state at two nearb...
thermo_t * m_thermo
pointer to the object representing the phase
DenseMatrix m_cstar
Dense matrix for cstar.
void update_T()
Update basic temperature-dependent quantities if the temperature has changed.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
Definition: GasTransport.h:339
const doublereal Pi
Pi.
Definition: ct_defs.h:51
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:556
doublereal m_temp
Current value of the temperature at which the properties in this object are calculated (Kelvin)...
Definition: GasTransport.h:301
std::vector< vector_fp > m_cstar_poly
Fit for cstar collision integral.
Definition: GasTransport.h:385
static const doublereal Min_C_Internal
Constant to compare dimensionless heat capacities against zero.
doublereal Frot(doublereal tr, doublereal sqtr)
The Parker temperature correction to the rotational collision number.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:97
int m_mode
Type of the polynomial fits to temperature.
Definition: GasTransport.h:256
virtual void init(ThermoPhase *thermo, int mode=0, int log_level=0)
Initialize a transport manager.
R poly8(D x, R *c)
Templated evaluation of a polynomial of order 8.
Definition: utilities.h:616
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:64
std::vector< vector_fp > m_omega22_poly
Fit for omega22 collision integral.
Definition: GasTransport.h:361
ThermoPhase object for the ideal gas equation of state - workhorse for Cantera (see Thermodynamic Pro...
R poly6(D x, R *c)
Templated evaluation of a polynomial of order 6.
Definition: utilities.h:604
vector_fp m_molefracs_last
Mole fraction vector from last L-matrix evaluation.
virtual void updateDiff_T()
Update the binary diffusion coefficients.
virtual void init(thermo_t *thermo, int mode=0, int log_level=0)
Initialize a transport manager.
vector_fp m_visc
vector of species viscosities (kg /m /s).
Definition: GasTransport.h:266
virtual void getMolarFluxes(const doublereal *const state1, const doublereal *const state2, const doublereal delta, doublereal *const fluxes)
Get the molar diffusional fluxes [kmol/m^2/s] of the species, given the thermodynamic state at two ne...
vector_fp m_spwork
work space length = m_kk
Definition: GasTransport.h:262
int solve(doublereal *b, size_t nrhs=1, size_t ldb=0)
Solves the Ax = b system returning x in the b spot.
void eval_L0000(const doublereal *const x)
Evaluate the L0000 matrices.
const doublereal * massFractions() const
Return a const pointer to the mass fraction array.
Definition: Phase.h:482
virtual void getMultiDiffCoeffs(const size_t ld, doublereal *const d)
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
std::vector< vector_fp > m_bstar_poly
Fit for bstar collision integral.
Definition: GasTransport.h:377
int factor()
Factors the A matrix, overwriting A.
std::vector< vector_fp > m_astar_poly
Fit for astar collision integral.
Definition: GasTransport.h:369
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
vector_fp m_eps
Lennard-Jones well-depth of the species in the current phase.
Definition: GasTransport.h:418
void updateThermal_T()
Update the temperature-dependent terms needed to compute the thermal conductivity and thermal diffusi...
std::vector< vector_int > m_poly
Indices for the (i,j) interaction in collision integral fits.
Definition: GasTransport.h:353
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:265
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:278
vector_fp m_zrot
Rotational relaxation number for each species.
Definition: GasTransport.h:391
doublereal temperature() const
Temperature (K).
Definition: Phase.h:602
DenseMatrix m_epsilon
The effective well depth for (i,j) collisions.
Definition: GasTransport.h:453
vector_fp m_crot
Dimensionless rotational heat capacity of each species.
Definition: GasTransport.h:399
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
doublereal m_sqrt_t
current value of temperature to 1/2 power
Definition: GasTransport.h:311
const doublereal Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:142
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:669
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual void getSpeciesFluxes(size_t ndim, const doublereal *const grad_T, size_t ldx, const doublereal *const grad_X, size_t ldf, doublereal *const fluxes)
Get the species diffusive mass fluxes wrt to the mass averaged velocity, given the gradients in mole ...
size_t m_nsp
Number of species.
doublereal m_logt
Current value of the log of the temperature.
Definition: GasTransport.h:314
bool m_abc_ok
Boolean indicating viscosity is up to date.
virtual void getCp_R_ref(doublereal *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
Definition: ThermoPhase.h:849
doublereal m_kbt
Current value of Boltzmann constant times the temperature (Joules)
Definition: GasTransport.h:304
void eval_L0010(const doublereal *const x)
Evaluate the L0010 matrices.
const doublereal Boltzmann
Boltzmann's constant [J/K].
Definition: ct_defs.h:76
DenseMatrix m_om22
Dense matrix for omega22.
vector_fp m_molefracs
Vector of species mole fractions.
Definition: GasTransport.h:237
Class GasTransport implements some functions and properties that are shared by the MixTransport and M...
Definition: GasTransport.h:19
void update_C()
Update basic concentration-dependent quantities if the concentrations have changed.
vector_fp m_mw
Local copy of the species molecular weights.
Definition: GasTransport.h:274