53 const char* Foam::chemkinReader::reactionTypeNames[4] =
57 "nonEquilibriumReversible",
61 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
65 "unimolecularFallOff",
66 "chemicallyActivatedBimolecular",
70 "unknownReactionRateType"
73 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
78 "unknownFallOffFunctionType"
81 void Foam::chemkinReader::initReactionKeywordTable()
83 reactionKeywordTable_.
insert(
"M", thirdBodyReactionType);
84 reactionKeywordTable_.
insert(
"LOW", unimolecularFallOffReactionType);
85 reactionKeywordTable_.
insert(
"HIGH", chemicallyActivatedBimolecularReactionType);
86 reactionKeywordTable_.
insert(
"TROE", TroeReactionType);
87 reactionKeywordTable_.
insert(
"SRI", SRIReactionType);
88 reactionKeywordTable_.
insert(
"LT", LandauTellerReactionType);
89 reactionKeywordTable_.
insert(
"RLT", reverseLandauTellerReactionType);
90 reactionKeywordTable_.
insert(
"JAN", JanevReactionType);
91 reactionKeywordTable_.
insert(
"FIT1", powerSeriesReactionRateType);
92 reactionKeywordTable_.
insert(
"HV", radiationActivatedReactionType);
93 reactionKeywordTable_.
insert(
"TDEP", speciesTempReactionType);
94 reactionKeywordTable_.
insert(
"EXCI", energyLossReactionType);
95 reactionKeywordTable_.
insert(
"MOME", plasmaMomentumTransfer);
96 reactionKeywordTable_.
insert(
"XSMI", collisionCrossSection);
97 reactionKeywordTable_.
insert(
"REV", nonEquilibriumReversibleReactionType);
98 reactionKeywordTable_.
insert(
"DUPLICATE", duplicateReactionType);
99 reactionKeywordTable_.
insert(
"DUP", duplicateReactionType);
100 reactionKeywordTable_.
insert(
"FORD", speciesOrderForward);
101 reactionKeywordTable_.
insert(
"RORD", speciesOrderReverse);
102 reactionKeywordTable_.
insert(
"UNITS", UnitsOfReaction);
103 reactionKeywordTable_.
insert(
"END", end);
107 Foam::scalar Foam::chemkinReader::molecularWeight
109 const List<specieElement>& specieComposition
114 forAll (specieComposition, i)
116 label nAtoms = specieComposition[i].nAtoms;
117 const word& elementName = specieComposition[i].elementName;
119 if (isotopeAtomicWts_.found(elementName))
121 molWt += nAtoms*isotopeAtomicWts_[elementName];
130 <<
"Unknown element " << elementName
131 <<
" on line " << lineNo_-1 <<
nl
132 <<
" specieComposition: " << specieComposition
141 void Foam::chemkinReader::checkCoeffs
144 const char* reactionRateName,
148 if (reactionCoeffs.size() != nCoeffs)
151 <<
"Wrong number of coefficients for the " << reactionRateName
152 <<
" rate expression on line "
153 << lineNo_-1 <<
", should be "
154 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl
155 <<
"Coefficients are "
156 << reactionCoeffs <<
nl
161 template<
class ReactionRateType>
162 void Foam::chemkinReader::addReactionType
164 const reactionType rType,
165 DynamicList<gasReaction::specieCoeffs>& lhs,
166 DynamicList<gasReaction::specieCoeffs>& rhs,
167 const ReactionRateType& rr
176 new IrreversibleReaction<gasThermoPhysics, ReactionRateType>
178 Reaction<gasThermoPhysics>
195 new ReversibleReaction<gasThermoPhysics, ReactionRateType>
197 Reaction<gasThermoPhysics>
215 <<
"Reaction type " << reactionTypeNames[rType]
216 <<
" on line " << lineNo_-1
217 <<
" not handled by this function"
223 <<
"Unknown reaction type " << rType
224 <<
" on line " << lineNo_-1
230 template<
template<
class,
class>
class PressureDependencyType>
231 void Foam::chemkinReader::addPressureDependentReaction
233 const reactionType rType,
234 const fallOffFunctionType fofType,
235 DynamicList<gasReaction::specieCoeffs>& lhs,
236 DynamicList<gasReaction::specieCoeffs>& rhs,
240 const HashTable<scalarList>& reactionCoeffsTable,
241 const scalar Afactor0,
242 const scalar AfactorInf,
246 checkCoeffs(k0Coeffs,
"k0", 3);
247 checkCoeffs(kInfCoeffs,
"kInf", 3);
257 PressureDependencyType
258 <ArrheniusReactionRate, LindemannFallOffFunction>
260 ArrheniusReactionRate
262 Afactor0*k0Coeffs[0],
266 ArrheniusReactionRate
268 AfactorInf*kInfCoeffs[0],
272 LindemannFallOffFunction(),
273 thirdBodyEfficiencies(speciesTable_, efficiencies)
283 reactionCoeffsTable[fallOffFunctionNames[fofType]]
286 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
288 FatalErrorIn(
"chemkinReader::addPressureDependentReaction")
289 <<
"Wrong number of coefficients for Troe rate expression"
290 " on line " << lineNo_-1 <<
", should be 3 or 4 but "
291 << TroeCoeffs.size() <<
" supplied." <<
nl
292 <<
"Coefficients are "
297 if (TroeCoeffs.size() == 3)
299 TroeCoeffs.setSize(4);
300 TroeCoeffs[3] = GREAT;
307 PressureDependencyType
308 <ArrheniusReactionRate, TroeFallOffFunction>
310 ArrheniusReactionRate
312 Afactor0*k0Coeffs[0],
316 ArrheniusReactionRate
318 AfactorInf*kInfCoeffs[0],
329 thirdBodyEfficiencies(speciesTable_, efficiencies)
339 reactionCoeffsTable[fallOffFunctionNames[fofType]]
342 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
344 FatalErrorIn(
"chemkinReader::addPressureDependentReaction")
345 <<
"Wrong number of coefficients for SRI rate expression"
346 " on line " << lineNo_-1 <<
", should be 3 or 5 but "
347 << SRICoeffs.size() <<
" supplied." <<
nl
348 <<
"Coefficients are "
353 if (SRICoeffs.size() == 3)
355 SRICoeffs.setSize(5);
364 PressureDependencyType
365 <ArrheniusReactionRate, SRIFallOffFunction>
367 ArrheniusReactionRate
369 Afactor0*k0Coeffs[0],
373 ArrheniusReactionRate
375 AfactorInf*kInfCoeffs[0],
387 thirdBodyEfficiencies(speciesTable_, efficiencies)
397 FatalErrorIn(
"chemkinReader::addPressureDependentReaction")
398 <<
"Fall-off function type "
399 << fallOffFunctionNames[fofType]
400 <<
" on line " << lineNo_-1
401 <<
" not implemented"
406 FatalErrorIn(
"chemkinReader::addPressureDependentReaction")
407 <<
"Unknown fall-off function type " << fofType
408 <<
" on line " << lineNo_-1
416 void Foam::chemkinReader::addReaction
418 DynamicList<gasReaction::specieCoeffs>& lhs,
419 DynamicList<gasReaction::specieCoeffs>& rhs,
421 const reactionType rType,
422 const reactionRateType rrType,
423 const fallOffFunctionType fofType,
425 HashTable<scalarList>& reactionCoeffsTable,
429 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
435 const List<specieElement>& specieComposition =
436 specieComposition_[speciesTable_[lhs[i].index]];
438 forAll (specieComposition, j)
440 label elementi = elementIndices_[specieComposition[j].elementName];
441 nAtoms[elementi] += lhs[i].stoichCoeff*specieComposition[j].nAtoms;
447 const List<specieElement>& specieComposition =
448 specieComposition_[speciesTable_[rhs[i].index]];
450 forAll (specieComposition, j)
452 label elementi = elementIndices_[specieComposition[j].elementName];
453 nAtoms[elementi] -= rhs[i].stoichCoeff*specieComposition[j].nAtoms;
460 const scalar concFactor = 0.001;
464 sumExp += lhs[i].exponent;
466 scalar Afactor =
pow(concFactor, sumExp - 1.0);
468 scalar AfactorRev = Afactor;
470 if (rType == nonEquilibriumReversible)
475 sumExp += rhs[i].exponent;
477 AfactorRev =
pow(concFactor, sumExp - 1.0);
484 if (rType == nonEquilibriumReversible)
487 reactionCoeffsTable[reactionTypeNames[rType]];
489 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
493 new NonEquilibriumReversibleReaction
496 Reaction<gasThermoPhysics>
503 ArrheniusReactionRate
505 Afactor*ArrheniusCoeffs[0],
507 ArrheniusCoeffs[2]/RR
509 ArrheniusReactionRate
511 AfactorRev*reverseArrheniusCoeffs[0],
512 reverseArrheniusCoeffs[1],
513 reverseArrheniusCoeffs[2]/RR
524 ArrheniusReactionRate
526 Afactor*ArrheniusCoeffs[0],
528 ArrheniusCoeffs[2]/RR
535 case thirdBodyArrhenius:
537 if (rType == nonEquilibriumReversible)
540 reactionCoeffsTable[reactionTypeNames[rType]];
542 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
546 new NonEquilibriumReversibleReaction
549 Reaction<gasThermoPhysics>
556 thirdBodyArrheniusReactionRate
558 Afactor*concFactor*ArrheniusCoeffs[0],
560 ArrheniusCoeffs[2]/RR,
561 thirdBodyEfficiencies(speciesTable_, efficiencies)
563 thirdBodyArrheniusReactionRate
565 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
566 reverseArrheniusCoeffs[1],
567 reverseArrheniusCoeffs[2]/RR,
568 thirdBodyEfficiencies(speciesTable_, efficiencies)
579 thirdBodyArrheniusReactionRate
581 Afactor*concFactor*ArrheniusCoeffs[0],
583 ArrheniusCoeffs[2]/RR,
584 thirdBodyEfficiencies(speciesTable_, efficiencies)
591 case unimolecularFallOff:
593 addPressureDependentReaction<FallOffReactionRate>
600 reactionCoeffsTable[reactionRateTypeNames[rrType]],
610 case chemicallyActivatedBimolecular:
612 addPressureDependentReaction<ChemicallyActivatedReactionRate>
620 reactionCoeffsTable[reactionRateTypeNames[rrType]],
632 reactionCoeffsTable[reactionRateTypeNames[rrType]];
633 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
635 if (rType == nonEquilibriumReversible)
638 reactionCoeffsTable[reactionTypeNames[rType]];
639 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
644 word(reactionTypeNames[rType])
645 + reactionRateTypeNames[rrType]
647 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
651 new NonEquilibriumReversibleReaction
654 Reaction<gasThermoPhysics>
661 LandauTellerReactionRate
663 Afactor*ArrheniusCoeffs[0],
665 ArrheniusCoeffs[2]/RR,
666 LandauTellerCoeffs[0],
667 LandauTellerCoeffs[1]
669 LandauTellerReactionRate
671 AfactorRev*reverseArrheniusCoeffs[0],
672 reverseArrheniusCoeffs[1],
673 reverseArrheniusCoeffs[2]/RR,
674 reverseLandauTellerCoeffs[0],
675 reverseLandauTellerCoeffs[1]
686 LandauTellerReactionRate
688 Afactor*ArrheniusCoeffs[0],
690 ArrheniusCoeffs[2]/RR,
691 LandauTellerCoeffs[0],
692 LandauTellerCoeffs[1]
702 reactionCoeffsTable[reactionRateTypeNames[rrType]];
704 checkCoeffs(JanevCoeffs,
"Janev", 9);
712 Afactor*ArrheniusCoeffs[0],
714 ArrheniusCoeffs[2]/RR,
724 reactionCoeffsTable[reactionRateTypeNames[rrType]];
726 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
732 powerSeriesReactionRate
734 Afactor*ArrheniusCoeffs[0],
736 ArrheniusCoeffs[2]/RR,
737 powerSeriesCoeffs.begin()
743 case unknownReactionRateType:
746 <<
"Internal error on line " << lineNo_-1
747 <<
": reaction rate type has not been set"
757 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
758 <<
" on line " << lineNo_-1
759 <<
" not implemented"
765 <<
"Unknown reaction rate type " << rrType
766 <<
" on line " << lineNo_-1
775 if (
mag(nAtoms[i]) > SMALL)
778 <<
"Elemental imbalance in " << elementNames_[i]
779 <<
" in reaction" <<
nl
780 << reactions_.last() <<
nl
781 <<
" on line " << lineNo_-1
788 reactionCoeffsTable.clear();
792 void Foam::chemkinReader::read
794 const fileName& CHEMKINFileName,
795 const fileName& thermoFileName
800 std::ifstream thermoStream(thermoFileName.c_str());
806 "chemkin::chemkin(const fileName& CHEMKINFileName, "
807 "const fileName& thermoFileName)"
808 ) <<
"file " << thermoFileName <<
" not found"
812 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
813 yy_switch_to_buffer(bufferPtr);
818 yy_delete_buffer(bufferPtr);
823 std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
829 "chemkin::chemkin(const fileName& CHEMKINFileName, "
830 "const fileName& thermoFileName)"
831 ) <<
"file " << CHEMKINFileName <<
" not found"
835 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
836 yy_switch_to_buffer(bufferPtr);
838 initReactionKeywordTable();
843 yy_delete_buffer(bufferPtr);
850 Foam::chemkinReader::chemkinReader
858 speciesTable_(static_cast<const wordList&>(
wordList()))
860 read(CHEMKINFileName, thermoFileName);
865 Foam::chemkinReader::chemkinReader(
const dictionary& thermoDict)
878 if (thermoDict.
found(
"CHEMKINThermoFile"))
880 thermoFile =
fileName(thermoDict.
lookup(
"CHEMKINThermoFile")).expand();
887 if (chemkinFile.size() && chemkinFile[0] !=
'/')
889 chemkinFile = relPath/chemkinFile;
892 if (thermoFile.size() && thermoFile[0] !=
'/')
894 thermoFile = relPath/thermoFile;
898 read(chemkinFile, thermoFile);