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