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