27 double Frot(
double tr,
double sqtr)
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;
37 void MultiTransport::init(
ThermoPhase* thermo,
int mode,
int log_level)
39 GasTransport::init(thermo, mode, log_level);
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);
56 m_lmatrix_soln_ok =
false;
57 m_thermal_tlast = 0.0;
60 m_spwork1.resize(m_nsp);
61 m_spwork2.resize(m_nsp);
62 m_spwork3.resize(m_nsp);
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);
75 const double sq298 = sqrt(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);
84 double MultiTransport::thermalConductivity()
86 solveLMatrixEquation();
88 for (
size_t k = 0; k < 2*m_nsp; k++) {
89 sum += m_b[k + m_nsp] * m_a[k + m_nsp];
94 void MultiTransport::getThermalDiffCoeffs(
double*
const dt)
96 solveLMatrixEquation();
98 for (
size_t k = 0; k < m_nsp; k++) {
99 dt[k] = c * m_mw[k] * m_molefracs[k] * m_a[k];
103 double MultiTransport::pressure_ig()
105 return m_thermo->molarDensity() *
GasConstant * m_thermo->temperature();
108 void MultiTransport::solveLMatrixEquation()
113 if (m_lmatrix_soln_ok) {
120 for (
size_t k = 0; k < m_nsp; k++) {
122 m_b[k + m_nsp] = m_molefracs[k];
123 m_b[k + 2*m_nsp] = m_molefracs[k];
135 for (
size_t k = 0; k < m_nsp; k++) {
136 if (!hasInternalModes(k)) {
137 m_b[2*m_nsp + k] = 0.0;
142 m_Lmatrix.resize(3*m_nsp, 3*m_nsp, 0.0);
145 eval_L0000(m_molefracs.data());
146 eval_L0010(m_molefracs.data());
149 eval_L1010(m_molefracs.data());
150 eval_L1001(m_molefracs.data());
153 eval_L0101(m_molefracs.data());
158 solve(m_Lmatrix, m_a.data());
159 m_lmatrix_soln_ok =
true;
160 m_molefracs_last = m_molefracs;
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)
175 bool addThermalDiffusion =
false;
176 for (
size_t i = 0; i < ndim; i++) {
177 if (grad_T[i] != 0.0) {
178 addThermalDiffusion =
true;
181 if (addThermalDiffusion) {
182 getThermalDiffCoeffs(m_spwork.data());
185 const double* y = m_thermo->massFractions();
186 double rho = m_thermo->density();
188 for (
size_t i = 0; i < m_nsp; i++) {
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);
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]);
211 for (
size_t j = 0; j < m_nsp; j++) {
214 vector<double> gsave(ndim), grx(ldx*m_nsp);
215 for (
size_t n = 0; n < ldx*ndim; n++) {
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;
227 solve(m_aa, fluxes, ndim, ldf);
228 double pp = pressure_ig();
232 for (
size_t n = 0; n < ndim; n++) {
234 for (
size_t i = 0; i < m_nsp; i++) {
235 fluxes[i +
offset] *= rho * y[i] / pp;
240 if (addThermalDiffusion) {
241 for (
size_t n = 0; n < ndim; n++) {
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;
251 void MultiTransport::getMassFluxes(
const double* state1,
const double* state2,
252 double delta,
double* fluxes)
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);
263 m_thermo->restoreState(nsp+2, state2);
264 double p2 = m_thermo->pressure();
265 double t2 = state2[0];
266 m_thermo->getMoleFractions(x2);
268 double p = 0.5*(p1 + p2);
269 double t = 0.5*(state1[0] + state2[0]);
271 for (
size_t n = 0; n < nsp; n++) {
272 x3[n] = 0.5*(x1[n] + x2[n]);
274 m_thermo->setState_TPX(t, p, x3);
275 m_thermo->getMoleFractions(m_molefracs.data());
283 bool addThermalDiffusion =
false;
284 if (state1[0] != state2[0]) {
285 addThermalDiffusion =
true;
286 getThermalDiffCoeffs(m_spwork.data());
289 const double* y = m_thermo->massFractions();
290 double rho = m_thermo->density();
291 for (
size_t i = 0; i < m_nsp; i++) {
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);
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]);
314 for (
size_t j = 0; j < m_nsp; j++) {
316 fluxes[j] = x2[j] - x1[j];
323 double pp = pressure_ig();
326 for (
size_t i = 0; i < m_nsp; i++) {
327 fluxes[i] *= rho * y[i] / pp;
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;
339 void MultiTransport::getMolarFluxes(
const double*
const state1,
340 const double*
const state2,
342 double*
const fluxes)
344 getMassFluxes(state1, state2, delta, fluxes);
345 for (
size_t k = 0; k < m_thermo->nSpecies(); k++) {
346 fluxes[k] /= m_mw[k];
350 void MultiTransport::getMultiDiffCoeffs(
const size_t ld,
double*
const d)
352 double p = pressure_ig();
364 eval_L0000(m_molefracs.data());
368 int ierr =
invert(m_Lmatrix, m_nsp);
370 throw CanteraError(
"MultiTransport::getMultiDiffCoeffs",
371 "invert returned ierr = {}", ierr);
374 m_lmatrix_soln_ok =
false;
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));
388 void MultiTransport::update_T()
390 if (m_temp == m_thermo->temperature() && m_nsp == m_thermo->nSpecies()) {
393 GasTransport::update_T();
396 m_lmatrix_soln_ok =
false;
400 void MultiTransport::update_C()
403 m_thermo->getMoleFractions(m_molefracs.data());
405 for (
size_t k = 0; k < m_nsp; k++) {
407 m_molefracs[k] = std::max(
Tiny, m_molefracs[k]);
408 if (m_molefracs[k] != m_molefracs_last[k]) {
412 m_lmatrix_soln_ok =
false;
417 void MultiTransport::updateThermal_T()
419 if (m_thermal_tlast == m_thermo->temperature()) {
423 updateSpeciesViscosities();
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());
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());
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);
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);
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];
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;
472 m_thermal_tlast = m_thermo->temperature();
478 bool MultiTransport::hasInternalModes(
size_t j)
483 void MultiTransport::eval_L0000(
const double*
const x)
485 double prefactor = 16.0*m_temp/25.0;
487 for (
size_t i = 0; i < m_nsp; i++) {
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);
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));
501 m_Lmatrix(i,i) = 0.0;
505 void MultiTransport::eval_L0010(
const double*
const x)
507 double prefactor = 1.6*m_temp;
508 for (
size_t j = 0; j < m_nsp; j++) {
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));
519 sum -= m_Lmatrix(i,j+m_nsp);
521 m_Lmatrix(j,j+m_nsp) += sum;
525 void MultiTransport::eval_L1000()
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);
534 void MultiTransport::eval_L1010(
const double* x)
536 const double fiveover3pi = 5.0/(3.0*
Pi);
537 double prefactor = (16.0*m_temp)/25.0;
539 for (
size_t j = 0; j < m_nsp; 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];
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])));
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)
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]);
564 m_Lmatrix(j+m_nsp,j+m_nsp) -= sum*constant1;
568 void MultiTransport::eval_L1001(
const double* x)
570 double prefactor = 32.00*m_temp/(5.00*
Pi);
571 for (
size_t j = 0; j < m_nsp; j++) {
573 if (hasInternalModes(j)) {
574 double constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]);
576 for (
size_t i = 0; i < m_nsp; i++) {
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);
582 m_Lmatrix(j+m_nsp,j+2*m_nsp) += sum;
584 for (
size_t i = 0; i < m_nsp; i++) {
585 m_Lmatrix(i+m_nsp,j+2*m_nsp) = 0.0;
591 void MultiTransport::eval_L0001()
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;
600 void MultiTransport::eval_L0100()
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;
609 void MultiTransport::eval_L0110()
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);
618 void MultiTransport::eval_L0101(
const double* x)
620 for (
size_t i = 0; i < m_nsp; i++) {
621 if (hasInternalModes(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]);
627 for (
size_t k = 0; k < m_nsp; k++) {
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;
633 sum += x[k]*m_astar(i,k)*constant2 / (m_mw[k]*diff_int);
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])
642 for (
size_t k = 0; k < m_nsp; k++) {
643 m_Lmatrix(i+2*m_nsp,i+2*m_nsp) = 1.0;
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.
Base class for a phase with thermodynamic properties.
R poly8(D x, R *c)
Templated evaluation of a polynomial of order 8.
R poly6(D x, R *c)
Templated evaluation of a polynomial of order 6.
const double Boltzmann
Boltzmann constant [J/K].
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
const double Tiny
Small number to compare differences of mole fractions against.
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.
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...