48 const char* Foam::chemkinReader::reactionTypeNames[4] =
52 "nonEquilibriumReversible",
56 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
60 "unimolecularFallOff",
61 "chemicallyActivatedBimolecular",
65 "unknownReactionRateType"
68 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
73 "unknownFallOffFunctionType"
76 void Foam::chemkinReader::initReactionKeywordTable()
78 reactionKeywordTable_.
insert(
"M", thirdBodyReactionType);
79 reactionKeywordTable_.
insert(
"LOW", unimolecularFallOffReactionType);
80 reactionKeywordTable_.
insert
83 chemicallyActivatedBimolecularReactionType
85 reactionKeywordTable_.
insert(
"TROE", TroeReactionType);
86 reactionKeywordTable_.
insert(
"SRI", SRIReactionType);
87 reactionKeywordTable_.
insert(
"LT", LandauTellerReactionType);
88 reactionKeywordTable_.
insert(
"RLT", reverseLandauTellerReactionType);
89 reactionKeywordTable_.
insert(
"JAN", JanevReactionType);
90 reactionKeywordTable_.
insert(
"FIT1", powerSeriesReactionRateType);
91 reactionKeywordTable_.
insert(
"HV", radiationActivatedReactionType);
92 reactionKeywordTable_.
insert(
"TDEP", speciesTempReactionType);
93 reactionKeywordTable_.
insert(
"EXCI", energyLossReactionType);
94 reactionKeywordTable_.
insert(
"MOME", plasmaMomentumTransfer);
95 reactionKeywordTable_.
insert(
"XSMI", collisionCrossSection);
96 reactionKeywordTable_.
insert(
"REV", nonEquilibriumReversibleReactionType);
97 reactionKeywordTable_.
insert(
"DUPLICATE", duplicateReactionType);
98 reactionKeywordTable_.
insert(
"DUP", duplicateReactionType);
99 reactionKeywordTable_.
insert(
"FORD", speciesOrderForward);
100 reactionKeywordTable_.
insert(
"RORD", speciesOrderReverse);
101 reactionKeywordTable_.
insert(
"UNITS", UnitsOfReaction);
102 reactionKeywordTable_.
insert(
"END", end);
106 Foam::scalar Foam::chemkinReader::molecularWeight
108 const List<specieElement>& specieComposition
113 forAll(specieComposition, i)
115 label nAtoms = specieComposition[i].nAtoms();
116 const word& elementName = specieComposition[i].name();
118 if (isotopeAtomicWts_.found(elementName))
120 molWt += nAtoms*isotopeAtomicWts_[elementName];
129 <<
"Unknown element " << elementName
130 <<
" on line " << lineNo_-1 <<
nl
131 <<
" specieComposition: " << specieComposition
140 void Foam::chemkinReader::checkCoeffs
143 const char* reactionRateName,
147 if (reactionCoeffs.size() != nCoeffs)
150 <<
"Wrong number of coefficients for the " << reactionRateName
151 <<
" rate expression on line "
152 << lineNo_-1 <<
", should be "
153 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl
154 <<
"Coefficients are "
155 << reactionCoeffs <<
nl
160 template<
class ReactionRateType>
161 void Foam::chemkinReader::addReactionType
163 const reactionType rType,
164 DynamicList<specieCoeffs>& lhs,
165 DynamicList<specieCoeffs>& rhs,
166 const ReactionRateType& rr
175 new IrreversibleReaction<thermoPhysics, ReactionRateType>
177 ReactionProxy<thermoPhysics>
194 new ReversibleReaction<thermoPhysics, ReactionRateType>
196 ReactionProxy<thermoPhysics>
214 <<
"Reaction type " << reactionTypeNames[rType]
215 <<
" on line " << lineNo_-1
216 <<
" not handled by this function"
222 <<
"Unknown reaction type " << rType
223 <<
" on line " << lineNo_-1
229 template<
template<
class,
class>
class PressureDependencyType>
230 void Foam::chemkinReader::addPressureDependentReaction
232 const reactionType rType,
233 const fallOffFunctionType fofType,
234 DynamicList<specieCoeffs>& lhs,
235 DynamicList<specieCoeffs>& rhs,
239 const HashTable<scalarList>& reactionCoeffsTable,
240 const scalar Afactor0,
241 const scalar AfactorInf,
245 checkCoeffs(k0Coeffs,
"k0", 3);
246 checkCoeffs(kInfCoeffs,
"kInf", 3);
256 PressureDependencyType
257 <ArrheniusReactionRate, LindemannFallOffFunction>
259 ArrheniusReactionRate
261 Afactor0*k0Coeffs[0],
265 ArrheniusReactionRate
267 AfactorInf*kInfCoeffs[0],
271 LindemannFallOffFunction(),
272 thirdBodyEfficiencies(speciesTable_, efficiencies)
281 reactionCoeffsTable[fallOffFunctionNames[fofType]]
284 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
287 <<
"Wrong number of coefficients for Troe rate expression"
288 " on line " << lineNo_-1 <<
", should be 3 or 4 but "
289 << TroeCoeffs.size() <<
" supplied." <<
nl
290 <<
"Coefficients are "
295 if (TroeCoeffs.size() == 3)
297 TroeCoeffs.setSize(4);
298 TroeCoeffs[3] = great;
305 PressureDependencyType
306 <ArrheniusReactionRate, TroeFallOffFunction>
308 ArrheniusReactionRate
310 Afactor0*k0Coeffs[0],
314 ArrheniusReactionRate
316 AfactorInf*kInfCoeffs[0],
327 thirdBodyEfficiencies(speciesTable_, efficiencies)
336 reactionCoeffsTable[fallOffFunctionNames[fofType]]
339 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
342 <<
"Wrong number of coefficients for SRI rate expression"
343 " on line " << lineNo_-1 <<
", should be 3 or 5 but "
344 << SRICoeffs.size() <<
" supplied." <<
nl
345 <<
"Coefficients are "
350 if (SRICoeffs.size() == 3)
352 SRICoeffs.setSize(5);
361 PressureDependencyType
362 <ArrheniusReactionRate, SRIFallOffFunction>
364 ArrheniusReactionRate
366 Afactor0*k0Coeffs[0],
370 ArrheniusReactionRate
372 AfactorInf*kInfCoeffs[0],
384 thirdBodyEfficiencies(speciesTable_, efficiencies)
392 <<
"Fall-off function type "
393 << fallOffFunctionNames[fofType]
394 <<
" on line " << lineNo_-1
395 <<
" not implemented"
402 void Foam::chemkinReader::addReaction
404 DynamicList<specieCoeffs>& lhs,
405 DynamicList<specieCoeffs>& rhs,
407 const reactionType rType,
408 const reactionRateType rrType,
409 const fallOffFunctionType fofType,
411 HashTable<scalarList>& reactionCoeffsTable,
415 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
421 const List<specieElement>& specieComposition =
422 speciesComposition_[speciesTable_[lhs[i].index]];
424 forAll(specieComposition, j)
426 label elementi = elementIndices_[specieComposition[j].name()];
428 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
434 const List<specieElement>& specieComposition =
435 speciesComposition_[speciesTable_[rhs[i].index]];
437 forAll(specieComposition, j)
439 label elementi = elementIndices_[specieComposition[j].name()];
441 rhs[i].stoichCoeff*specieComposition[j].nAtoms();
448 const scalar concFactor = 0.001;
452 sumExp += lhs[i].exponent;
454 scalar Afactor =
pow(concFactor, sumExp - 1.0);
456 scalar AfactorRev = Afactor;
458 if (rType == nonEquilibriumReversible)
463 sumExp += rhs[i].exponent;
465 AfactorRev =
pow(concFactor, sumExp - 1.0);
472 if (rType == nonEquilibriumReversible)
475 reactionCoeffsTable[reactionTypeNames[rType]];
477 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
481 new NonEquilibriumReversibleReaction
484 ArrheniusReactionRate
487 ReactionProxy<thermoPhysics>
494 ArrheniusReactionRate
496 Afactor*ArrheniusCoeffs[0],
498 ArrheniusCoeffs[2]/
RR
500 ArrheniusReactionRate
502 AfactorRev*reverseArrheniusCoeffs[0],
503 reverseArrheniusCoeffs[1],
504 reverseArrheniusCoeffs[2]/
RR
515 ArrheniusReactionRate
517 Afactor*ArrheniusCoeffs[0],
519 ArrheniusCoeffs[2]/
RR
525 case thirdBodyArrhenius:
527 if (rType == nonEquilibriumReversible)
530 reactionCoeffsTable[reactionTypeNames[rType]];
532 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
536 new NonEquilibriumReversibleReaction
539 thirdBodyArrheniusReactionRate
542 ReactionProxy<thermoPhysics>
549 thirdBodyArrheniusReactionRate
551 Afactor*concFactor*ArrheniusCoeffs[0],
553 ArrheniusCoeffs[2]/
RR,
554 thirdBodyEfficiencies(speciesTable_, efficiencies)
556 thirdBodyArrheniusReactionRate
558 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
559 reverseArrheniusCoeffs[1],
560 reverseArrheniusCoeffs[2]/
RR,
561 thirdBodyEfficiencies(speciesTable_, efficiencies)
572 thirdBodyArrheniusReactionRate
574 Afactor*concFactor*ArrheniusCoeffs[0],
576 ArrheniusCoeffs[2]/
RR,
577 thirdBodyEfficiencies(speciesTable_, efficiencies)
583 case unimolecularFallOff:
585 addPressureDependentReaction<FallOffReactionRate>
592 reactionCoeffsTable[reactionRateTypeNames[rrType]],
601 case chemicallyActivatedBimolecular:
603 addPressureDependentReaction<ChemicallyActivatedReactionRate>
611 reactionCoeffsTable[reactionRateTypeNames[rrType]],
622 reactionCoeffsTable[reactionRateTypeNames[rrType]];
623 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
625 if (rType == nonEquilibriumReversible)
628 reactionCoeffsTable[reactionTypeNames[rType]];
629 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
634 word(reactionTypeNames[rType])
635 + reactionRateTypeNames[rrType]
637 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
641 new NonEquilibriumReversibleReaction
644 LandauTellerReactionRate
647 ReactionProxy<thermoPhysics>
654 LandauTellerReactionRate
656 Afactor*ArrheniusCoeffs[0],
658 ArrheniusCoeffs[2]/
RR,
659 LandauTellerCoeffs[0],
660 LandauTellerCoeffs[1]
662 LandauTellerReactionRate
664 AfactorRev*reverseArrheniusCoeffs[0],
665 reverseArrheniusCoeffs[1],
666 reverseArrheniusCoeffs[2]/
RR,
667 reverseLandauTellerCoeffs[0],
668 reverseLandauTellerCoeffs[1]
679 LandauTellerReactionRate
681 Afactor*ArrheniusCoeffs[0],
683 ArrheniusCoeffs[2]/
RR,
684 LandauTellerCoeffs[0],
685 LandauTellerCoeffs[1]
694 reactionCoeffsTable[reactionRateTypeNames[rrType]];
696 checkCoeffs(JanevCoeffs,
"Janev", 9);
704 Afactor*ArrheniusCoeffs[0],
706 ArrheniusCoeffs[2]/
RR,
707 FixedList<scalar, 9>(JanevCoeffs)
715 reactionCoeffsTable[reactionRateTypeNames[rrType]];
717 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
723 powerSeriesReactionRate
725 Afactor*ArrheniusCoeffs[0],
727 FixedList<scalar, 4>(powerSeriesCoeffs)
732 case unknownReactionRateType:
735 <<
"Internal error on line " << lineNo_-1
736 <<
": reaction rate type has not been set"
743 <<
"Reaction rate type index " << rrType
744 <<
" on line " << lineNo_-1
753 if (
mag(nAtoms[i]) > imbalanceTol_)
756 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
757 <<
" in " << elementNames_[i]
758 <<
" in reaction" <<
nl
759 << reactions_.last() <<
nl
760 <<
" on line " << lineNo_-1
767 reactionCoeffsTable.clear();
771 void Foam::chemkinReader::read
773 const fileName& CHEMKINFileName,
774 const fileName& thermoFileName,
775 const fileName& transportFileName
781 transportDict_.read(IFstream(transportFileName)());
785 std::ifstream thermoStream(thermoFileName.c_str());
790 <<
"file " << thermoFileName <<
" not found"
794 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
795 yy_switch_to_buffer(bufferPtr);
800 yy_delete_buffer(bufferPtr);
805 std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
810 <<
"file " << CHEMKINFileName <<
" not found"
814 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
815 yy_switch_to_buffer(bufferPtr);
817 initReactionKeywordTable();
822 yy_delete_buffer(bufferPtr);
828 Foam::chemkinReader::chemkinReader
830 const fileName& CHEMKINFileName,
831 const fileName& transportFileName,
832 const fileName& thermoFileName,
839 newFormat_(newFormat),
840 imbalanceTol_(rootSmall)
842 read(CHEMKINFileName, thermoFileName, transportFileName);
#define forAll(list, i)
Loop across all elements in list.
Macros for easy insertion into run-time selection tables.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
bool found(const Key &) const
Return true if hashedEntry is found in table.
static scalar TlowDefault
Default temperature limits of applicability of reaction rates.
static scalar ThighDefault
static const fileName null
An empty fileName.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
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.
List< scalar > scalarList
A List of scalars.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
atomicWeightTable atomicWeights