26 doublereal
Frot(doublereal tr, doublereal sqtr)
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;
36 MultiTransport::MultiTransport(
thermo_t* thermo)
47 m_a.resize(3*
m_nsp, 1.0);
48 m_b.resize(3*
m_nsp, 0.0);
51 m_frot_298.resize(
m_nsp);
52 m_rotrelax.resize(
m_nsp);
53 m_cinternal.resize(
m_nsp);
62 m_lmatrix_soln_ok =
false;
63 m_thermal_tlast = 0.0;
66 m_spwork1.resize(
m_nsp);
67 m_spwork2.resize(
m_nsp);
68 m_spwork3.resize(
m_nsp);
72 for (
size_t i = 0; i <
m_nsp; i++) {
73 for (
size_t j = i; j <
m_nsp; j++) {
75 m_log_eps_k(j,i) = m_log_eps_k(i,j);
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++) {
86 m_frot_298[k] =
Frot(
m_eps[k]/kb298, m_sqrt_eps_k[k]/sq298);
92 solveLMatrixEquation();
94 for (
size_t k = 0; k < 2*
m_nsp; k++) {
102 solveLMatrixEquation();
104 for (
size_t k = 0; k <
m_nsp; k++) {
109 void MultiTransport::solveLMatrixEquation()
114 if (m_lmatrix_soln_ok) {
121 for (
size_t k = 0; k <
m_nsp; k++) {
136 for (
size_t k = 0; k <
m_nsp; k++) {
137 if (!hasInternalModes(k)) {
138 m_b[2*
m_nsp + k] = 0.0;
159 solve(m_Lmatrix, m_a.data());
160 m_lmatrix_soln_ok =
true;
167 size_t ldx,
const doublereal*
const grad_X,
168 size_t ldf, doublereal*
const fluxes)
176 bool addThermalDiffusion =
false;
177 for (
size_t i = 0; i < ndim; i++) {
178 if (grad_T[i] != 0.0) {
179 addThermalDiffusion =
true;
182 if (addThermalDiffusion) {
189 for (
size_t i = 0; i <
m_nsp; i++) {
191 for (
size_t j = 0; j <
m_nsp; j++) {
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]);
212 for (
size_t j = 0; j <
m_nsp; j++) {
216 for (
size_t n = 0; n < ldx*ndim; n++) {
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;
228 solve(m_aa, fluxes, ndim, ldf);
229 doublereal pp = pressure_ig();
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;
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;
255 double* x1 = m_spwork1.data();
256 double* x2 = m_spwork2.data();
257 double* x3 = m_spwork3.data();
261 double t1 = state1[0];
266 double t2 = state2[0];
269 double p = 0.5*(p1 + p2);
270 double t = 0.5*(state1[0] + state2[0]);
272 for (
size_t n = 0; n < nsp; n++) {
273 x3[n] = 0.5*(x1[n] + x2[n]);
284 bool addThermalDiffusion =
false;
285 if (state1[0] != state2[0]) {
286 addThermalDiffusion =
true;
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++) {
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]);
315 for (
size_t j = 0; j <
m_nsp; j++) {
317 fluxes[j] = x2[j] - x1[j];
324 doublereal pp = pressure_ig();
327 for (
size_t i = 0; i <
m_nsp; i++) {
328 fluxes[i] *= rho * y[i] / pp;
332 if (addThermalDiffusion) {
333 doublereal grad_logt = (t2 - t1)/
m_temp;
334 for (
size_t i = 0; i <
m_nsp; i++) {
341 const doublereal*
const state2,
342 const doublereal delta,
343 doublereal*
const fluxes)
347 fluxes[k] /=
m_mw[k];
353 doublereal p = pressure_ig();
371 throw CanteraError(
"MultiTransport::getMultiDiffCoeffs",
372 "invert returned ierr = {}", ierr);
375 m_lmatrix_soln_ok =
false;
377 doublereal prefactor = 16.0 *
m_temp 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];
383 (m_Lmatrix(i,j) - m_Lmatrix(i,i));
393 GasTransport::update_T();
397 m_lmatrix_soln_ok =
false;
406 for (
size_t k = 0; k <
m_nsp; k++) {
413 m_lmatrix_soln_ok =
false;
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);
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;
483 bool MultiTransport::hasInternalModes(
size_t j)
490 doublereal 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 doublereal 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++) {
539 void MultiTransport::eval_L1010(
const doublereal* x)
541 const doublereal fiveover3pi = 5.0/(3.0*
Pi);
542 doublereal 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) *
573 void MultiTransport::eval_L1001(
const doublereal* x)
575 doublereal 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++) {
596 void 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;
605 void 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;
614 void MultiTransport::eval_L0110()
616 for (
size_t j = 0; j <
m_nsp; j++) {
617 for (
size_t i = 0; i <
m_nsp; i++) {
623 void MultiTransport::eval_L0101(
const doublereal* 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++) {
virtual doublereal thermalConductivity()
Returns the mixture thermal conductivity in W/m/K.
void eval_L1000()
Evaluate the L1000 matrices.
virtual void updateSpeciesViscosities()
Update the pure-species viscosities.
DenseMatrix m_bstar
Dense matrix for bstar.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
void restoreState(const vector_fp &state)
Restore a state saved on a previous call to saveState.
virtual void getThermalDiffCoeffs(doublereal *const dt)
Return the thermal diffusion coefficients (kg/m/s)
doublereal temperature() const
Temperature (K).
const doublereal * massFractions() const
Return a const pointer to the mass fraction array.
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
Interface for class MultiTransport.
DenseMatrix m_astar
Dense matrix for astar.
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...
thermo_t * m_thermo
pointer to the object representing the phase
DenseMatrix m_cstar
Dense matrix for cstar.
size_t nSpecies() const
Returns the number of species in the phase.
virtual doublereal density() const
Density (kg/m^3).
void update_T()
Update basic temperature-dependent quantities if the temperature has changed.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
doublereal m_temp
Current value of the temperature at which the properties in this object are calculated (Kelvin)...
std::vector< vector_fp > m_cstar_poly
Fit for cstar collision integral.
static const doublereal Min_C_Internal
Constant to compare dimensionless heat capacities against zero.
virtual void getCp_R_ref(doublereal *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
doublereal Frot(doublereal tr, doublereal sqtr)
The Parker temperature correction to the rotational collision number.
Base class for a phase with thermodynamic properties.
int m_mode
Type of the polynomial fits to temperature.
virtual void init(ThermoPhase *thermo, int mode=0, int log_level=0)
Initialize a transport manager.
R poly8(D x, R *c)
Templated evaluation of a polynomial of order 8.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
std::vector< vector_fp > m_omega22_poly
Fit for omega22 collision integral.
ThermoPhase object for the ideal gas equation of state - workhorse for Cantera (see Thermodynamic Pro...
R poly6(D x, R *c)
Templated evaluation of a polynomial of order 6.
vector_fp m_molefracs_last
Mole fraction vector from last L-matrix evaluation.
virtual void updateDiff_T()
Update the binary diffusion coefficients.
virtual void init(thermo_t *thermo, int mode=0, int log_level=0)
Initialize a transport manager.
vector_fp m_visc
vector of species viscosities (kg /m /s).
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_spwork
work space length = m_kk
void eval_L0000(const doublereal *const x)
Evaluate the L0000 matrices.
virtual void getMultiDiffCoeffs(const size_t ld, doublereal *const d)
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
std::vector< vector_fp > m_bstar_poly
Fit for bstar collision integral.
std::vector< vector_fp > m_astar_poly
Fit for astar collision integral.
Base class for exceptions thrown by Cantera classes.
vector_fp m_eps
Lennard-Jones well-depth of the species in the current phase.
void updateThermal_T()
Update the temperature-dependent terms needed to compute the thermal conductivity and thermal diffusi...
std::vector< vector_int > m_poly
Indices for the (i,j) interaction in collision integral fits.
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
vector_fp m_zrot
Rotational relaxation number for each species.
DenseMatrix m_epsilon
The effective well depth for (i,j) collisions.
vector_fp m_crot
Dimensionless rotational heat capacity of each species.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
doublereal m_sqrt_t
current value of temperature to 1/2 power
const doublereal Tiny
Small number to compare differences of mole fractions against.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Contains declarations for string manipulation functions within Cantera.
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 ...
size_t m_nsp
Number of species.
doublereal m_logt
Current value of the log of the temperature.
bool m_abc_ok
Boolean indicating viscosity is up to date.
doublereal m_kbt
Current value of Boltzmann constant times the temperature (Joules)
void eval_L0010(const doublereal *const x)
Evaluate the L0010 matrices.
Namespace for the Cantera kernel.
const doublereal Boltzmann
Boltzmann's constant [J/K].
DenseMatrix m_om22
Dense matrix for omega22.
vector_fp m_molefracs
Vector of species mole fractions.
Class GasTransport implements some functions and properties that are shared by the MixTransport and M...
void update_C()
Update basic concentration-dependent quantities if the concentrations have changed.
vector_fp m_mw
Local copy of the species molecular weights.