27double 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;
48 m_a.resize(3*
m_nsp, 1.0);
49 m_b.resize(3*
m_nsp, 0.0);
52 m_frot_298.resize(
m_nsp);
53 m_rotrelax.resize(
m_nsp);
54 m_cinternal.resize(
m_nsp);
61 m_lmatrix_soln_ok =
false;
62 m_thermal_tlast = 0.0;
65 m_spwork1.resize(
m_nsp);
66 m_spwork2.resize(
m_nsp);
67 m_spwork3.resize(
m_nsp);
71 for (
size_t i = 0; i <
m_nsp; i++) {
72 for (
size_t j = i; j <
m_nsp; j++) {
74 m_log_eps_k(j,i) = m_log_eps_k(i,j);
80 const double sq298 = sqrt(298.0);
82 m_sqrt_eps_k.resize(
m_nsp);
83 for (
size_t k = 0; k <
m_nsp; k++) {
85 m_frot_298[k] =
Frot(
m_eps[k]/kb298, m_sqrt_eps_k[k]/sq298);
91 solveLMatrixEquation();
93 for (
size_t k = 0; k < 2*
m_nsp; k++) {
101 solveLMatrixEquation();
103 for (
size_t k = 0; k <
m_nsp; k++) {
108double MultiTransport::pressure_ig()
113void MultiTransport::solveLMatrixEquation()
118 if (m_lmatrix_soln_ok) {
125 for (
size_t k = 0; k <
m_nsp; k++) {
140 for (
size_t k = 0; k <
m_nsp; k++) {
141 if (!hasInternalModes(k)) {
142 m_b[2*
m_nsp + k] = 0.0;
163 solve(m_Lmatrix, m_a.data());
164 m_lmatrix_soln_ok =
true;
171 size_t ldx,
const double*
const grad_X,
172 size_t ldf,
double*
const fluxes)
180 bool addThermalDiffusion =
false;
181 for (
size_t i = 0; i < ndim; i++) {
182 if (grad_T[i] != 0.0) {
183 addThermalDiffusion =
true;
186 if (addThermalDiffusion) {
193 for (
size_t i = 0; i <
m_nsp; i++) {
195 for (
size_t j = 0; j <
m_nsp; j++) {
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]);
216 for (
size_t j = 0; j <
m_nsp; j++) {
219 vector<double> gsave(ndim), grx(ldx*
m_nsp);
220 for (
size_t n = 0; n < ldx*ndim; n++) {
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;
232 solve(m_aa, fluxes, ndim, ldf);
233 double pp = pressure_ig();
237 for (
size_t n = 0; n < ndim; n++) {
239 for (
size_t i = 0; i <
m_nsp; i++) {
240 fluxes[i +
offset] *= rho * y[i] / pp;
245 if (addThermalDiffusion) {
246 for (
size_t n = 0; n < ndim; n++) {
248 double grad_logt = grad_T[n]/
m_temp;
249 for (
size_t i = 0; i <
m_nsp; i++) {
257 double delta,
double* fluxes)
259 double* x1 = m_spwork1.data();
260 double* x2 = m_spwork2.data();
261 double* x3 = m_spwork3.data();
265 double t1 = state1[0];
270 double t2 = state2[0];
273 double p = 0.5*(p1 + p2);
274 double t = 0.5*(state1[0] + state2[0]);
276 for (
size_t n = 0; n < nsp; n++) {
277 x3[n] = 0.5*(x1[n] + x2[n]);
288 bool addThermalDiffusion =
false;
289 if (state1[0] != state2[0]) {
290 addThermalDiffusion =
true;
296 for (
size_t i = 0; i <
m_nsp; i++) {
298 for (
size_t j = 0; j <
m_nsp; j++) {
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]);
319 for (
size_t j = 0; j <
m_nsp; j++) {
321 fluxes[j] = x2[j] - x1[j];
328 double pp = pressure_ig();
331 for (
size_t i = 0; i <
m_nsp; i++) {
332 fluxes[i] *= rho * y[i] / pp;
336 if (addThermalDiffusion) {
337 double grad_logt = (t2 - t1)/
m_temp;
338 for (
size_t i = 0; i <
m_nsp; i++) {
345 const double*
const state2,
347 double*
const fluxes)
351 fluxes[k] /=
m_mw[k];
357 double p = pressure_ig();
375 throw CanteraError(
"MultiTransport::getMultiDiffCoeffs",
376 "invert returned ierr = {}", ierr);
379 m_lmatrix_soln_ok =
false;
381 double prefactor = 16.0 *
m_temp
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];
387 (m_Lmatrix(i,j) - m_Lmatrix(i,i));
398 GasTransport::update_T();
401 m_lmatrix_soln_ok =
false;
410 for (
size_t k = 0; k <
m_nsp; k++) {
417 m_lmatrix_soln_ok =
false;
432 for (
size_t i = 0; i <
m_nsp; i++) {
433 for (
size_t j = i; j <
m_nsp; j++) {
452 for (
size_t k = 0; k <
m_nsp; k++) {
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);
459 for (
size_t k = 0; k <
m_nsp; k++) {
474 for (
size_t k = 0; k <
m_nsp; k++) {
475 m_cinternal[k] = cp[k] - 2.5;
483bool MultiTransport::hasInternalModes(
size_t j)
490 double prefactor = 16.0*
m_temp/25.0;
492 for (
size_t i = 0; i <
m_nsp; i++) {
496 for (
size_t k = 0; k <
m_nsp; k++) {
501 for (
size_t j = 0; j !=
m_nsp; ++j) {
502 m_Lmatrix(i,j) = prefactor * x[j]
506 m_Lmatrix(i,i) = 0.0;
512 double prefactor = 1.6*
m_temp;
513 for (
size_t j = 0; j <
m_nsp; j++) {
517 for (
size_t i = 0; i <
m_nsp; i++) {
518 m_Lmatrix(i,j +
m_nsp) = - prefactor * x[i] * xj *
m_mw[i] *
524 sum -= m_Lmatrix(i,j+
m_nsp);
526 m_Lmatrix(j,j+
m_nsp) += sum;
532 for (
size_t j = 0; j <
m_nsp; j++) {
533 for (
size_t i = 0; i <
m_nsp; i++) {
539void MultiTransport::eval_L1010(
const double* x)
541 const double fiveover3pi = 5.0/(3.0*
Pi);
542 double prefactor = (16.0*
m_temp)/25.0;
544 for (
size_t j = 0; j <
m_nsp; j++) {
546 double constant1 = prefactor*x[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];
554 for (
size_t i = 0; i <
m_nsp; i++) {
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])));
561 (constant2 - threemjsq*
m_bstar(i,j)
564 sum += x[i] /(term1) *
573void MultiTransport::eval_L1001(
const double* x)
575 double prefactor = 32.00*
m_temp/(5.00*
Pi);
576 for (
size_t j = 0; j <
m_nsp; j++) {
578 if (hasInternalModes(j)) {
579 double constant = prefactor*
m_mw[j]*x[j]*
m_crot[j]/(m_cinternal[j]*m_rotrelax[j]);
581 for (
size_t i = 0; i <
m_nsp; i++) {
589 for (
size_t i = 0; i <
m_nsp; i++) {
596void MultiTransport::eval_L0001()
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;
605void MultiTransport::eval_L0100()
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;
614void MultiTransport::eval_L0110()
616 for (
size_t j = 0; j <
m_nsp; j++) {
617 for (
size_t i = 0; i <
m_nsp; i++) {
623void MultiTransport::eval_L0101(
const double* x)
625 for (
size_t i = 0; i <
m_nsp; i++) {
626 if (hasInternalModes(i)) {
628 double constant1 = 4*
m_temp*x[i]/m_cinternal[i];
630 (5*
Pi*m_cinternal[i]*m_rotrelax[i]);
632 for (
size_t k = 0; k <
m_nsp; k++) {
634 double diff_int =
m_bdiff(i,k);
636 sum += x[k]/diff_int;
638 sum += x[k]*
m_astar(i,k)*constant2 / (
m_mw[k]*diff_int);
647 for (
size_t k = 0; k <
m_nsp; k++) {
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).
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
size_t nSpecies() const
Returns the number of species in the phase.
double temperature() const
Temperature (K).
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
const double * massFractions() const
Return a const pointer to the mass fraction array.
virtual double density() const
Density (kg/m^3).
virtual double pressure() const
Return the thermodynamic pressure (Pa).
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
size_t m_nsp
Number of species.
ThermoPhase & thermo()
Phase object.
R poly6(D x, R *c)
Templated evaluation of a polynomial of order 6.
R poly8(D x, R *c)
Templated evaluation of a polynomial of order 8.
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.
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...