Cantera  3.0.0
Loading...
Searching...
No Matches
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
13
14using namespace std;
15
16namespace 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 */
27double 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
38 : GasTransport(thermo)
39{
40}
41
42void 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 m_frot_298.resize(m_nsp);
53 m_rotrelax.resize(m_nsp);
54 m_cinternal.resize(m_nsp);
58
59 // set flags all false
60 m_l0000_ok = false;
61 m_lmatrix_soln_ok = false;
62 m_thermal_tlast = 0.0;
63
64 // some work space
65 m_spwork1.resize(m_nsp);
66 m_spwork2.resize(m_nsp);
67 m_spwork3.resize(m_nsp);
68
69 // precompute and store log(epsilon_ij/k_B)
70 m_log_eps_k.resize(m_nsp, m_nsp);
71 for (size_t i = 0; i < m_nsp; i++) {
72 for (size_t j = i; j < m_nsp; j++) {
73 m_log_eps_k(i,j) = log(m_epsilon(i,j)/Boltzmann);
74 m_log_eps_k(j,i) = m_log_eps_k(i,j);
75 }
76 }
77
78 // precompute and store constant parts of the Parker rotational
79 // collision number temperature correction
80 const double sq298 = sqrt(298.0);
81 const double kb298 = Boltzmann * 298.0;
82 m_sqrt_eps_k.resize(m_nsp);
83 for (size_t k = 0; k < m_nsp; k++) {
84 m_sqrt_eps_k[k] = sqrt(m_eps[k]/Boltzmann);
85 m_frot_298[k] = Frot(m_eps[k]/kb298, m_sqrt_eps_k[k]/sq298);
86 }
87}
88
90{
91 solveLMatrixEquation();
92 double sum = 0.0;
93 for (size_t k = 0; k < 2*m_nsp; k++) {
94 sum += m_b[k + m_nsp] * m_a[k + m_nsp];
95 }
96 return -4.0*sum;
97}
98
100{
101 solveLMatrixEquation();
102 const double c = 1.6/GasConstant;
103 for (size_t k = 0; k < m_nsp; k++) {
104 dt[k] = c * m_mw[k] * m_molefracs[k] * m_a[k];
105 }
106}
107
108double MultiTransport::pressure_ig()
109{
111}
112
113void MultiTransport::solveLMatrixEquation()
114{
115 // if T has changed, update the temperature-dependent properties.
117 update_C();
118 if (m_lmatrix_soln_ok) {
119 return;
120 }
121
122 // Copy the mole fractions twice into the last two blocks of the right-hand-
123 // side vector m_b. The first block of m_b was set to zero when it was
124 // created, and is not modified so doesn't need to be reset to zero.
125 for (size_t k = 0; k < m_nsp; k++) {
126 m_b[k] = 0.0;
127 m_b[k + m_nsp] = m_molefracs[k];
128 m_b[k + 2*m_nsp] = m_molefracs[k];
129 }
130
131 // Set the right-hand side vector to zero in the 3rd block for all species
132 // with no internal energy modes. The corresponding third-block rows and
133 // columns will be set to zero, except on the diagonal of L01,01, where they
134 // are set to 1.0. This has the effect of eliminating these equations from
135 // the system, since the equation becomes: m_a[2*m_nsp + k] = 0.0.
136
137 // Note that this differs from the Chemkin procedure, where all *monatomic*
138 // species are excluded. Since monatomic radicals can have non-zero internal
139 // heat capacities due to electronic excitation, they should be retained.
140 for (size_t k = 0; k < m_nsp; k++) {
141 if (!hasInternalModes(k)) {
142 m_b[2*m_nsp + k] = 0.0;
143 }
144 }
145
146 // evaluate the submatrices of the L matrix
147 m_Lmatrix.resize(3*m_nsp, 3*m_nsp, 0.0);
148
149 //! Evaluate the upper-left block of the L matrix.
150 eval_L0000(m_molefracs.data());
151 eval_L0010(m_molefracs.data());
152 eval_L0001();
153 eval_L1000();
154 eval_L1010(m_molefracs.data());
155 eval_L1001(m_molefracs.data());
156 eval_L0100();
157 eval_L0110();
158 eval_L0101(m_molefracs.data());
159
160 // Solve it using GMRES or LU decomposition. The last solution in m_a should
161 // provide a good starting guess, so convergence should be fast.
162 m_a = m_b;
163 solve(m_Lmatrix, m_a.data());
164 m_lmatrix_soln_ok = true;
166 // L matrix is overwritten with LU decomposition
167 m_l0000_ok = false;
168}
169
170void MultiTransport::getSpeciesFluxes(size_t ndim, const double* const grad_T,
171 size_t ldx, const double* const grad_X,
172 size_t ldf, double* const fluxes)
173{
174 // update the binary diffusion coefficients if necessary
175 update_T();
176 updateDiff_T();
177
178 // If any component of grad_T is non-zero, then get the
179 // thermal diffusion coefficients
180 bool addThermalDiffusion = false;
181 for (size_t i = 0; i < ndim; i++) {
182 if (grad_T[i] != 0.0) {
183 addThermalDiffusion = true;
184 }
185 }
186 if (addThermalDiffusion) {
188 }
189
190 const double* y = m_thermo->massFractions();
191 double rho = m_thermo->density();
192
193 for (size_t i = 0; i < m_nsp; i++) {
194 double sum = 0.0;
195 for (size_t j = 0; j < m_nsp; j++) {
196 m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j);
197 sum += m_aa(i,j);
198 }
199 m_aa(i,i) -= sum;
200 }
201
202 // enforce the condition \sum Y_k V_k = 0. This is done by replacing
203 // the flux equation with the largest gradx component in the first
204 // coordinate direction with the flux balance condition.
205 size_t jmax = 0;
206 double gradmax = -1.0;
207 for (size_t j = 0; j < m_nsp; j++) {
208 if (fabs(grad_X[j]) > gradmax) {
209 gradmax = fabs(grad_X[j]);
210 jmax = j;
211 }
212 }
213
214 // set the matrix elements in this row to the mass fractions,
215 // and set the entry in gradx to zero
216 for (size_t j = 0; j < m_nsp; j++) {
217 m_aa(jmax,j) = y[j];
218 }
219 vector<double> gsave(ndim), grx(ldx*m_nsp);
220 for (size_t n = 0; n < ldx*ndim; n++) {
221 grx[n] = grad_X[n];
222 }
223
224 // copy grad_X to fluxes
225 for (size_t n = 0; n < ndim; n++) {
226 const double* gx = grad_X + ldx*n;
227 copy(gx, gx + m_nsp, fluxes + ldf*n);
228 fluxes[jmax + n*ldf] = 0.0;
229 }
230
231 // solve the equations
232 solve(m_aa, fluxes, ndim, ldf);
233 double pp = pressure_ig();
234
235 // multiply diffusion velocities by rho * V to create mass fluxes, and
236 // restore the gradx elements that were modified
237 for (size_t n = 0; n < ndim; n++) {
238 size_t offset = n*ldf;
239 for (size_t i = 0; i < m_nsp; i++) {
240 fluxes[i + offset] *= rho * y[i] / pp;
241 }
242 }
243
244 // thermal diffusion
245 if (addThermalDiffusion) {
246 for (size_t n = 0; n < ndim; n++) {
247 size_t offset = n*ldf;
248 double grad_logt = grad_T[n]/m_temp;
249 for (size_t i = 0; i < m_nsp; i++) {
250 fluxes[i + offset] -= m_spwork[i]*grad_logt;
251 }
252 }
253 }
254}
255
256void MultiTransport::getMassFluxes(const double* state1, const double* state2,
257 double delta, double* fluxes)
258{
259 double* x1 = m_spwork1.data();
260 double* x2 = m_spwork2.data();
261 double* x3 = m_spwork3.data();
262 size_t nsp = m_thermo->nSpecies();
263 m_thermo->restoreState(nsp+2, state1);
264 double p1 = m_thermo->pressure();
265 double t1 = state1[0];
267
268 m_thermo->restoreState(nsp+2, state2);
269 double p2 = m_thermo->pressure();
270 double t2 = state2[0];
272
273 double p = 0.5*(p1 + p2);
274 double t = 0.5*(state1[0] + state2[0]);
275
276 for (size_t n = 0; n < nsp; n++) {
277 x3[n] = 0.5*(x1[n] + x2[n]);
278 }
279 m_thermo->setState_TPX(t, p, x3);
281
282 // update the binary diffusion coefficients if necessary
283 update_T();
284 updateDiff_T();
285
286 // If there is a temperature gradient, then get the
287 // thermal diffusion coefficients
288 bool addThermalDiffusion = false;
289 if (state1[0] != state2[0]) {
290 addThermalDiffusion = true;
292 }
293
294 const double* y = m_thermo->massFractions();
295 double rho = m_thermo->density();
296 for (size_t i = 0; i < m_nsp; i++) {
297 double sum = 0.0;
298 for (size_t j = 0; j < m_nsp; j++) {
299 m_aa(i,j) = m_molefracs[j]*m_molefracs[i]/m_bdiff(i,j);
300 sum += m_aa(i,j);
301 }
302 m_aa(i,i) -= sum;
303 }
304
305 // enforce the condition \sum Y_k V_k = 0. This is done by replacing the
306 // flux equation with the largest gradx component with the flux balance
307 // condition.
308 size_t jmax = 0;
309 double gradmax = -1.0;
310 for (size_t j = 0; j < m_nsp; j++) {
311 if (fabs(x2[j] - x1[j]) > gradmax) {
312 gradmax = fabs(x1[j] - x2[j]);
313 jmax = j;
314 }
315 }
316
317 // set the matrix elements in this row to the mass fractions,
318 // and set the entry in gradx to zero
319 for (size_t j = 0; j < m_nsp; j++) {
320 m_aa(jmax,j) = y[j];
321 fluxes[j] = x2[j] - x1[j];
322 }
323 fluxes[jmax] = 0.0;
324
325 // Solve the equations
326 solve(m_aa, fluxes);
327
328 double pp = pressure_ig();
329 // multiply diffusion velocities by rho * Y_k to create
330 // mass fluxes, and divide by pressure
331 for (size_t i = 0; i < m_nsp; i++) {
332 fluxes[i] *= rho * y[i] / pp;
333 }
334
335 // thermal diffusion
336 if (addThermalDiffusion) {
337 double grad_logt = (t2 - t1)/m_temp;
338 for (size_t i = 0; i < m_nsp; i++) {
339 fluxes[i] -= m_spwork[i]*grad_logt;
340 }
341 }
342}
343
344void MultiTransport::getMolarFluxes(const double* const state1,
345 const double* const state2,
346 const double delta,
347 double* const fluxes)
348{
349 getMassFluxes(state1, state2, delta, fluxes);
350 for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
351 fluxes[k] /= m_mw[k];
352 }
353}
354
355void MultiTransport::getMultiDiffCoeffs(const size_t ld, double* const d)
356{
357 double p = pressure_ig();
358
359 // update the mole fractions
360 update_C();
361
362 // update the binary diffusion coefficients
363 update_T();
365
366 // evaluate L0000 if the temperature or concentrations have
367 // changed since it was last evaluated.
368 if (!m_l0000_ok) {
369 eval_L0000(m_molefracs.data());
370 }
371
372 // invert L00,00
373 int ierr = invert(m_Lmatrix, m_nsp);
374 if (ierr != 0) {
375 throw CanteraError("MultiTransport::getMultiDiffCoeffs",
376 "invert returned ierr = {}", ierr);
377 }
378 m_l0000_ok = false; // matrix is overwritten by inverse
379 m_lmatrix_soln_ok = false;
380
381 double prefactor = 16.0 * m_temp
382 * m_thermo->meanMolecularWeight()/(25.0 * p);
383 for (size_t i = 0; i < m_nsp; i++) {
384 for (size_t j = 0; j < m_nsp; j++) {
385 double c = prefactor/m_mw[j];
386 d[ld*j + i] = c*m_molefracs[i]*
387 (m_Lmatrix(i,j) - m_Lmatrix(i,i));
388 }
389 }
390
391}
392
394{
395 if (m_temp == m_thermo->temperature() && m_nsp == m_thermo->nSpecies()) {
396 return;
397 }
398 GasTransport::update_T();
399 // temperature has changed, so polynomial fits will need to be
400 // redone, and the L matrix reevaluated.
401 m_lmatrix_soln_ok = false;
402 m_l0000_ok = false;
403}
404
406{
407 // Update the local mole fraction array
409
410 for (size_t k = 0; k < m_nsp; k++) {
411 // add an offset to avoid a pure species condition
412 m_molefracs[k] = std::max(Tiny, m_molefracs[k]);
413 if (m_molefracs[k] != m_molefracs_last[k]) {
414 // If any mole fractions have changed, signal that concentration-
415 // dependent quantities will need to be recomputed before use.
416 m_l0000_ok = false;
417 m_lmatrix_soln_ok = false;
418 }
419 }
420}
421
423{
424 if (m_thermal_tlast == m_thermo->temperature()) {
425 return;
426 }
427 // we need species viscosities and binary diffusion coefficients
429 updateDiff_T();
430
431 // evaluate polynomial fits for A*, B*, C*
432 for (size_t i = 0; i < m_nsp; i++) {
433 for (size_t j = i; j < m_nsp; j++) {
434 double z = m_star_poly_uses_actualT[i][j] == 1 ? m_logt : m_logt - m_log_eps_k(i,j);
435 int ipoly = m_poly[i][j];
436 if (m_mode == CK_Mode) {
437 m_astar(i,j) = poly6(z, m_astar_poly[ipoly].data());
438 m_bstar(i,j) = poly6(z, m_bstar_poly[ipoly].data());
439 m_cstar(i,j) = poly6(z, m_cstar_poly[ipoly].data());
440 } else {
441 m_astar(i,j) = poly8(z, m_astar_poly[ipoly].data());
442 m_bstar(i,j) = poly8(z, m_bstar_poly[ipoly].data());
443 m_cstar(i,j) = poly8(z, m_cstar_poly[ipoly].data());
444 }
445 m_astar(j,i) = m_astar(i,j);
446 m_bstar(j,i) = m_bstar(i,j);
447 m_cstar(j,i) = m_cstar(i,j);
448 }
449 }
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 double 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<double> 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
481static const double Min_C_Internal = 0.001;
482
483bool MultiTransport::hasInternalModes(size_t j)
484{
485 return (m_cinternal[j] > Min_C_Internal);
486}
487
488void MultiTransport::eval_L0000(const double* const x)
489{
490 double prefactor = 16.0*m_temp/25.0;
491 double 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
510void MultiTransport::eval_L0010(const double* const x)
511{
512 double 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
539void MultiTransport::eval_L1010(const double* x)
540{
541 const double fiveover3pi = 5.0/(3.0*Pi);
542 double 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
573void MultiTransport::eval_L1001(const double* x)
574{
575 double 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
596void 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
605void 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
614void 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
623void MultiTransport::eval_L0101(const double* 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.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
Class GasTransport implements some functions and properties that are shared by the MixTransport and M...
vector< double > m_mw
Local copy of the species molecular weights.
vector< double > m_molefracs
Vector of species mole fractions.
double m_temp
Current value of the temperature at which the properties in this object are calculated (Kelvin).
vector< double > m_eps
Lennard-Jones well-depth of the species in the current phase.
virtual void updateDiff_T()
Update the binary diffusion coefficients.
DenseMatrix m_epsilon
The effective well depth for (i,j) collisions.
vector< double > m_spwork
work space length = m_kk
vector< double > m_zrot
Rotational relaxation number for each species.
int m_mode
Type of the polynomial fits to temperature.
double m_logt
Current value of the log of the temperature.
vector< double > m_visc
vector of species viscosities (kg /m /s).
virtual void updateSpeciesViscosities()
Update the pure-species viscosities.
double m_sqrt_t
current value of temperature to 1/2 power
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
vector< double > m_crot
Dimensionless rotational heat capacity of each species.
double m_kbt
Current value of Boltzmann constant times the temperature (Joules)
vector< vector< int > > m_star_poly_uses_actualT
Flag to indicate for which (i,j) interaction pairs the actual temperature is used instead of the redu...
vector< vector< double > > m_bstar_poly
Fit for bstar collision integral.
vector< vector< double > > m_astar_poly
Fit for astar collision integral.
vector< vector< int > > m_poly
Indices for the (i,j) interaction in collision integral fits.
vector< vector< double > > m_cstar_poly
Fit for cstar collision integral.
void init(ThermoPhase *thermo, int mode=0, int log_level=0) override
Initialize a transport manager.
void getMassFluxes(const double *state1, const double *state2, double delta, double *fluxes) override
Get the mass diffusional fluxes [kg/m^2/s] of the species, given the thermodynamic state at two nearb...
MultiTransport(ThermoPhase *thermo=0)
default constructor
void getMolarFluxes(const double *const state1, const double *const state2, const double delta, double *const fluxes) override
Get the molar diffusional fluxes [kmol/m^2/s] of the species, given the thermodynamic state at two ne...
bool m_l0000_ok
Boolean indicating viscosity is up to date.
void update_T() override
Update basic temperature-dependent quantities if the temperature has changed.
double thermalConductivity() override
Returns the mixture thermal conductivity in W/m/K.
void eval_L0010(const double *const x)
Evaluate the L0010 matrices.
void eval_L0000(const double *const x)
Evaluate the L0000 matrices.
DenseMatrix m_astar
Dense matrix for astar.
void getSpeciesFluxes(size_t ndim, const double *const grad_T, size_t ldx, const double *const grad_X, size_t ldf, double *const fluxes) override
Get the species diffusive mass fluxes wrt to the mass averaged velocity, given the gradients in mole ...
void updateThermal_T()
Update the temperature-dependent terms needed to compute the thermal conductivity and thermal diffusi...
vector< double > m_molefracs_last
Mole fraction vector from last L-matrix evaluation.
DenseMatrix m_cstar
Dense matrix for cstar.
void update_C() override
Update basic concentration-dependent quantities if the concentrations have changed.
void getThermalDiffCoeffs(double *const dt) override
Return the thermal diffusion coefficients (kg/m/s)
void getMultiDiffCoeffs(const size_t ld, double *const d) override
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
void init(ThermoPhase *thermo, int mode=0, int log_level=0) override
Initialize a transport manager.
DenseMatrix m_bstar
Dense matrix for bstar.
void eval_L1000()
Evaluate the L1000 matrices.
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:689
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:275
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
double temperature() const
Temperature (K).
Definition Phase.h:662
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:760
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:540
const double * massFractions() const
Return a const pointer to the mass fraction array.
Definition Phase.h:535
virtual double density() const
Density (kg/m^3).
Definition Phase.h:687
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:680
Base class for a phase with thermodynamic properties.
virtual void getCp_R_ref(double *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
virtual void setState_TPX(double t, double p, const double *x)
Set the temperature (K), pressure (Pa), and mole fractions.
ThermoPhase * m_thermo
pointer to the object representing the phase
Definition Transport.h:821
size_t m_nsp
Number of species.
Definition Transport.h:828
ThermoPhase & thermo()
Phase object.
Definition Transport.h:192
R poly6(D x, R *c)
Templated evaluation of a polynomial of order 6.
Definition utilities.h:117
R poly8(D x, R *c)
Templated evaluation of a polynomial of order 8.
Definition utilities.h:129
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.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
double Frot(double tr, double sqtr)
The Parker temperature correction to the rotational collision number.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...