31 template<
class ThermoType>
34 template<
class ThermoType>
40 template<
class ThermoType>
43 const PtrList<ThermoType>& speciesThermo
46 typename ThermoType::thermoType rhsThermo
49 *speciesThermo[rhs()[0].index].
W()
50 *speciesThermo[rhs()[0].index]
53 for (
label i=1; i<rhs().size(); ++i)
57 *speciesThermo[rhs()[i].index].W()
58 *speciesThermo[rhs()[i].index];
61 typename ThermoType::thermoType lhsThermo
64 *speciesThermo[lhs()[0].index].
W()
65 *speciesThermo[lhs()[0].index]
68 for (
label i=1; i<lhs().size(); ++i)
72 *speciesThermo[lhs()[i].index].W()
73 *speciesThermo[lhs()[i].index];
78 if (
mag(lhsThermo.Y() - rhsThermo.Y()) > 0.1)
81 <<
"Mass imbalance for reaction " <<
name() <<
": "
82 <<
mag(lhsThermo.Y() - rhsThermo.Y()) <<
" kg/kmol"
86 ThermoType::thermoType::operator=(lhsThermo == rhsThermo);
92 template<
class ThermoType>
102 ThermoType::thermoType(speciesThermo[0]),
106 setThermo(speciesThermo);
110 template<
class ThermoType>
118 ThermoType::thermoType(r),
124 template<
class ThermoType>
133 ThermoType::thermoType(speciesThermo[0]),
134 Tlow_(
dict.lookupOrDefault<scalar>(
"Tlow", TlowDefault)),
135 Thigh_(
dict.lookupOrDefault<scalar>(
"Thigh", ThighDefault))
137 setThermo(speciesThermo);
143 template<
class ThermoType>
154 typename dictionaryConstructorTable::iterator cstrIter =
155 dictionaryConstructorTablePtr_->find(reactionTypeName);
162 if (cstrIter == dictionaryConstructorTablePtr_->end())
164 cstrIter = dictionaryConstructorTablePtr_->find
170 if (cstrIter == dictionaryConstructorTablePtr_->end())
173 <<
"Unknown reaction type "
174 << reactionTypeName <<
nl <<
nl
175 <<
"Valid reaction types are :" <<
nl
176 << dictionaryConstructorTablePtr_->sortedToc()
182 cstrIter()(species, speciesThermo,
dict)
187 template<
class ThermoType>
199 if (!objectRegistryConstructorTablePtr_)
201 return New(species, speciesThermo,
dict);
206 typename objectRegistryConstructorTable::iterator cstrIter =
207 objectRegistryConstructorTablePtr_->find(reactionTypeName);
210 if (cstrIter == objectRegistryConstructorTablePtr_->end())
212 cstrIter = objectRegistryConstructorTablePtr_->find
218 if (cstrIter == objectRegistryConstructorTablePtr_->end())
220 typename dictionaryConstructorTable::iterator cstrIter =
221 dictionaryConstructorTablePtr_->find(reactionTypeName);
224 if (cstrIter == dictionaryConstructorTablePtr_->end())
226 cstrIter = dictionaryConstructorTablePtr_->find
232 if (cstrIter == dictionaryConstructorTablePtr_->end())
235 <<
"Unknown reaction type "
236 << reactionTypeName <<
nl <<
nl
237 <<
"Valid reaction types are :" <<
nl
238 << dictionaryConstructorTablePtr_->sortedToc()
239 << objectRegistryConstructorTablePtr_->sortedToc()
245 cstrIter()(species, speciesThermo,
dict)
251 cstrIter()(species, speciesThermo, ob,
dict)
258 template<
class ThermoType>
265 template<
class ThermoType>
280 const label si = lhs()[i].index;
282 Cf *=
c[si] >= small || el >= 1 ?
pow(
max(
c[si], 0), el) : 0;
287 const label si = rhs()[i].index;
289 Cr *=
c[si] >= small || er >= 1 ?
pow(
max(
c[si], 0), er) : 0;
294 template<
class ThermoType>
305 const scalar clippedT =
min(
max(
T, this->Tlow()), this->Thigh());
308 const scalar kf = this->kf(
p, clippedT,
c, li);
309 const scalar kr = this->kr(kf,
p, clippedT,
c, li);
313 this->
C(p,
T,
c, li, Cf, Cr);
318 return omegaf - omegar;
322 template<
class ThermoType>
335 scalar omegaf, omegar;
336 const scalar
omega = this->
omega(p,
T,
c, li, omegaf, omegar);
340 const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
341 const scalar sl = lhs()[i].stoichCoeff;
342 dNdtByV[Nsi0 + si] -= sl*
omega;
346 const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
347 const scalar sr = rhs()[i].stoichCoeff;
348 dNdtByV[Nsi0 + si] += sr*
omega;
353 template<
class ThermoType>
371 const scalar kf = this->kf(
p,
T,
c, li);
372 const scalar kr = this->kr(kf,
p,
T,
c, li);
376 this->
C(p,
T,
c, li, Cf, Cr);
379 const scalar
omega = kf*Cf - kr*Cr;
384 const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
385 const scalar sl = lhs()[i].stoichCoeff;
386 dNdtByV[Nsi0 + si] -= sl*
omega;
390 const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
391 const scalar sr = rhs()[i].stoichCoeff;
392 dNdtByV[Nsi0 + si] += sr*
omega;
400 const label sj = reduced ? c2s[lhs()[j].index] : lhs()[j].index;
405 const label si = lhs()[i].index;
410 c[si] >= small || el >= 1
417 c[si] >= small || el >= 1
425 const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
426 const scalar sl = lhs()[i].stoichCoeff;
427 ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sl*kf*dCfdcj;
431 const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
432 const scalar sr = rhs()[i].stoichCoeff;
433 ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sr*kf*dCfdcj;
439 const label sj = reduced ? c2s[rhs()[j].index] : rhs()[j].index;
444 const label si = rhs()[i].index;
449 c[si] >= small || er >= 1
456 c[si] >= small || er >= 1
464 const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
465 const scalar sl = lhs()[i].stoichCoeff;
466 ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sl*kr*dCrcj;
470 const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
471 const scalar sr = rhs()[i].stoichCoeff;
472 ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sr*kr*dCrcj;
480 const scalar dkfdT = this->dkfdT(
p,
T,
c, li);
481 const scalar dkrdT = this->dkrdT(
p,
T,
c, li, dkfdT, kr);
483 const scalar dwdT = dkfdT*Cf - dkrdT*Cr;
486 const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
487 const scalar sl = lhs()[i].stoichCoeff;
488 ddNdtByVdcTp(Nsi0 + si, Tsi) -= sl*dwdT;
492 const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
493 const scalar sr = rhs()[i].stoichCoeff;
494 ddNdtByVdcTp(Nsi0 + si, Tsi) += sr*dwdT;
505 this->dkfdc(
p,
T,
c, li, dkfdc);
506 this->dkrdc(
p,
T,
c, li, dkfdc, kr, dkrdc);
510 const label sj = reduced ? c2s[j] : j;
512 if (sj == -1)
continue;
514 const scalar dwdc = dkfdc[j]*Cf - dkrdc[j]*Cr;
517 const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
518 const scalar sl = lhs()[i].stoichCoeff;
519 ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sl*dwdc;
523 const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
524 const scalar sr = rhs()[i].stoichCoeff;
525 ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sr*dwdc;
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
#define forAll(list, i)
Loop across all elements in list.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Simple extension of ThermoType to handle reaction kinetics in addition to the equilibrium thermodynam...
virtual void write(Ostream &) const
Write.
static scalar TlowDefault
Default temperature limits of applicability of reaction rates.
static scalar ThighDefault
void ddNdtByVdcTp(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dNdtByV, scalarSquareMatrix &ddNdtByVdcTp, const bool reduced, const List< label > &c2s, const label csi0, const label Tsi, scalarField &cTpWork0, scalarField &cTpWork1) const
Derivative of the net reaction rate for each species involved.
void C(const scalar p, const scalar T, const scalarField &c, const label li, scalar &Cf, scalar &Cr) const
Concentration powers.
void dNdtByV(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dNdtByV, const bool reduced, const List< label > &c2s, const label Nsi0) const
The net reaction rate for each species involved.
static autoPtr< Reaction< ThermoType > > New(const speciesTable &species, const PtrList< ThermoType > &speciesThermo, const dictionary &dict)
Return a pointer to new reaction created from a dictionary.
scalar omega(const scalar p, const scalar T, const scalarField &c, const label li, scalar &omegaf, scalar &omegar) const
Net reaction rate.
Reaction(const speciesTable &species, const PtrList< ThermoType > &speciesThermo, const List< specieCoeffs > &lhs, const List< specieCoeffs > &rhs)
Construct from components.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
A list of keyword definitions, which are a keyword followed by any number of values (e....
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
A wordList with hashed indices for faster lookup by name.
Registry of regIOobjects.
Reaction base-class holding the specie names and coefficients.
void write(Ostream &) const
Write.
Type for exponents of species in reaction. Can take either integer or scalar value,...
bool removeTrailing(const char)
Remove trailing character returning true if string changed.
A class for handling words, derived from string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const dimensionedScalar c
Speed of light in a vacuum.
errorManipArg< error, int > exit(error &err, const int errNo=1)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
word name(const complex &)
Return a string representation of a complex.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const scalarList W(::W(thermo))