58 const char* Foam::chemkinReader::reactionTypeNames[4] =
62 "nonEquilibriumReversible",
66 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
70 "unimolecularFallOff",
71 "chemicallyActivatedBimolecular",
75 "unknownReactionRateType"
78 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
83 "unknownFallOffFunctionType"
86 void Foam::chemkinReader::initReactionKeywordTable()
88 reactionKeywordTable_.
insert(
"M", thirdBodyReactionType);
89 reactionKeywordTable_.
insert(
"LOW", unimolecularFallOffReactionType);
90 reactionKeywordTable_.
insert
93 chemicallyActivatedBimolecularReactionType
95 reactionKeywordTable_.
insert(
"TROE", TroeReactionType);
96 reactionKeywordTable_.
insert(
"SRI", SRIReactionType);
97 reactionKeywordTable_.
insert(
"LT", LandauTellerReactionType);
98 reactionKeywordTable_.
insert(
"RLT", reverseLandauTellerReactionType);
99 reactionKeywordTable_.
insert(
"JAN", JanevReactionType);
100 reactionKeywordTable_.
insert(
"FIT1", powerSeriesReactionRateType);
101 reactionKeywordTable_.
insert(
"HV", radiationActivatedReactionType);
102 reactionKeywordTable_.
insert(
"TDEP", speciesTempReactionType);
103 reactionKeywordTable_.
insert(
"EXCI", energyLossReactionType);
104 reactionKeywordTable_.
insert(
"MOME", plasmaMomentumTransfer);
105 reactionKeywordTable_.
insert(
"XSMI", collisionCrossSection);
106 reactionKeywordTable_.
insert(
"REV", nonEquilibriumReversibleReactionType);
107 reactionKeywordTable_.
insert(
"DUPLICATE", duplicateReactionType);
108 reactionKeywordTable_.
insert(
"DUP", duplicateReactionType);
109 reactionKeywordTable_.
insert(
"FORD", speciesOrderForward);
110 reactionKeywordTable_.
insert(
"RORD", speciesOrderReverse);
111 reactionKeywordTable_.
insert(
"UNITS", UnitsOfReaction);
112 reactionKeywordTable_.
insert(
"END", end);
116 Foam::scalar Foam::chemkinReader::molecularWeight
118 const List<specieElement>& specieComposition
123 forAll(specieComposition, i)
125 label nAtoms = specieComposition[i].nAtoms();
126 const word& elementName = specieComposition[i].name();
128 if (isotopeAtomicWts_.found(elementName))
130 molWt += nAtoms*isotopeAtomicWts_[elementName];
139 <<
"Unknown element " << elementName
140 <<
" on line " << lineNo_-1 <<
nl
141 <<
" specieComposition: " << specieComposition
150 void Foam::chemkinReader::checkCoeffs
153 const char* reactionRateName,
157 if (reactionCoeffs.size() != nCoeffs)
160 <<
"Wrong number of coefficients for the " << reactionRateName
161 <<
" rate expression on line "
162 << lineNo_-1 <<
", should be "
163 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl
164 <<
"Coefficients are "
165 << reactionCoeffs <<
nl
170 template<
class ReactionRateType>
171 void Foam::chemkinReader::addReactionType
173 const reactionType rType,
174 DynamicList<gasHReaction::specieCoeffs>& lhs,
175 DynamicList<gasHReaction::specieCoeffs>& rhs,
176 const ReactionRateType& rr
185 new IrreversibleReaction
188 Reaction<gasHThermoPhysics>
205 new ReversibleReaction
208 Reaction<gasHThermoPhysics>
226 <<
"Reaction type " << reactionTypeNames[rType]
227 <<
" on line " << lineNo_-1
228 <<
" not handled by this function"
234 <<
"Unknown reaction type " << rType
235 <<
" on line " << lineNo_-1
241 template<
template<
class,
class>
class PressureDependencyType>
242 void Foam::chemkinReader::addPressureDependentReaction
244 const reactionType rType,
245 const fallOffFunctionType fofType,
246 DynamicList<gasHReaction::specieCoeffs>& lhs,
247 DynamicList<gasHReaction::specieCoeffs>& rhs,
251 const HashTable<scalarList>& reactionCoeffsTable,
252 const scalar Afactor0,
253 const scalar AfactorInf,
257 checkCoeffs(k0Coeffs,
"k0", 3);
258 checkCoeffs(kInfCoeffs,
"kInf", 3);
268 PressureDependencyType
269 <ArrheniusReactionRate, LindemannFallOffFunction>
271 ArrheniusReactionRate
273 Afactor0*k0Coeffs[0],
277 ArrheniusReactionRate
279 AfactorInf*kInfCoeffs[0],
283 LindemannFallOffFunction(),
284 thirdBodyEfficiencies(speciesTable_, efficiencies)
293 reactionCoeffsTable[fallOffFunctionNames[fofType]]
296 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
299 <<
"Wrong number of coefficients for Troe rate expression"
300 " on line " << lineNo_-1 <<
", should be 3 or 4 but "
301 << TroeCoeffs.size() <<
" supplied." <<
nl
302 <<
"Coefficients are "
307 if (TroeCoeffs.size() == 3)
309 TroeCoeffs.setSize(4);
310 TroeCoeffs[3] = GREAT;
317 PressureDependencyType
318 <ArrheniusReactionRate, TroeFallOffFunction>
320 ArrheniusReactionRate
322 Afactor0*k0Coeffs[0],
326 ArrheniusReactionRate
328 AfactorInf*kInfCoeffs[0],
339 thirdBodyEfficiencies(speciesTable_, efficiencies)
348 reactionCoeffsTable[fallOffFunctionNames[fofType]]
351 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
354 <<
"Wrong number of coefficients for SRI rate expression"
355 " on line " << lineNo_-1 <<
", should be 3 or 5 but "
356 << SRICoeffs.size() <<
" supplied." <<
nl
357 <<
"Coefficients are "
362 if (SRICoeffs.size() == 3)
364 SRICoeffs.setSize(5);
373 PressureDependencyType
374 <ArrheniusReactionRate, SRIFallOffFunction>
376 ArrheniusReactionRate
378 Afactor0*k0Coeffs[0],
382 ArrheniusReactionRate
384 AfactorInf*kInfCoeffs[0],
396 thirdBodyEfficiencies(speciesTable_, efficiencies)
404 <<
"Fall-off function type "
405 << fallOffFunctionNames[fofType]
406 <<
" on line " << lineNo_-1
407 <<
" not implemented"
414 void Foam::chemkinReader::addReaction
416 DynamicList<gasHReaction::specieCoeffs>& lhs,
417 DynamicList<gasHReaction::specieCoeffs>& rhs,
419 const reactionType rType,
420 const reactionRateType rrType,
421 const fallOffFunctionType fofType,
423 HashTable<scalarList>& reactionCoeffsTable,
427 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
433 const List<specieElement>& specieComposition =
434 speciesComposition_[speciesTable_[lhs[i].index]];
436 forAll(specieComposition, j)
438 label elementi = elementIndices_[specieComposition[j].name()];
440 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
446 const List<specieElement>& specieComposition =
447 speciesComposition_[speciesTable_[rhs[i].index]];
449 forAll(specieComposition, j)
451 label elementi = elementIndices_[specieComposition[j].name()];
453 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<gasHThermoPhysics>
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
534 case thirdBodyArrhenius:
536 if (rType == nonEquilibriumReversible)
539 reactionCoeffsTable[reactionTypeNames[rType]];
541 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
545 new NonEquilibriumReversibleReaction
549 thirdBodyArrheniusReactionRate
552 Reaction<gasHThermoPhysics>
559 thirdBodyArrheniusReactionRate
561 Afactor*concFactor*ArrheniusCoeffs[0],
563 ArrheniusCoeffs[2]/
RR,
564 thirdBodyEfficiencies(speciesTable_, efficiencies)
566 thirdBodyArrheniusReactionRate
568 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
569 reverseArrheniusCoeffs[1],
570 reverseArrheniusCoeffs[2]/
RR,
571 thirdBodyEfficiencies(speciesTable_, efficiencies)
582 thirdBodyArrheniusReactionRate
584 Afactor*concFactor*ArrheniusCoeffs[0],
586 ArrheniusCoeffs[2]/
RR,
587 thirdBodyEfficiencies(speciesTable_, efficiencies)
593 case unimolecularFallOff:
595 addPressureDependentReaction<FallOffReactionRate>
602 reactionCoeffsTable[reactionRateTypeNames[rrType]],
611 case chemicallyActivatedBimolecular:
613 addPressureDependentReaction<ChemicallyActivatedReactionRate>
621 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<gasHThermoPhysics>
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]
701 reactionCoeffsTable[reactionRateTypeNames[rrType]];
703 checkCoeffs(JanevCoeffs,
"Janev", 9);
711 Afactor*ArrheniusCoeffs[0],
713 ArrheniusCoeffs[2]/
RR,
714 FixedList<scalar, 9>(JanevCoeffs)
722 reactionCoeffsTable[reactionRateTypeNames[rrType]];
724 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
730 powerSeriesReactionRate
732 Afactor*ArrheniusCoeffs[0],
734 ArrheniusCoeffs[2]/
RR,
735 FixedList<scalar, 4>(powerSeriesCoeffs)
740 case unknownReactionRateType:
743 <<
"Internal error on line " << lineNo_-1
744 <<
": reaction rate type has not been set"
751 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
752 <<
" on line " << lineNo_-1
753 <<
" not implemented"
761 if (
mag(nAtoms[i]) > imbalanceTol_)
764 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
765 <<
" in " << elementNames_[i]
766 <<
" in reaction" <<
nl
767 << reactions_.last() <<
nl
768 <<
" on line " << lineNo_-1
775 reactionCoeffsTable.clear();
779 void Foam::chemkinReader::read
781 const fileName& CHEMKINFileName,
782 const fileName& thermoFileName,
783 const fileName& transportFileName
786 transportDict_.read(IFstream(transportFileName)());
788 if (!thermoFileName.empty())
790 std::ifstream thermoStream(thermoFileName);
795 <<
"file " << thermoFileName <<
" not found"
799 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
800 yy_switch_to_buffer(bufferPtr);
805 yy_delete_buffer(bufferPtr);
810 std::ifstream CHEMKINStream(CHEMKINFileName);
815 <<
"file " << CHEMKINFileName <<
" not found"
819 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
820 yy_switch_to_buffer(bufferPtr);
822 initReactionKeywordTable();
827 yy_delete_buffer(bufferPtr);
833 Foam::chemkinReader::chemkinReader
844 speciesTable_(species),
845 reactions_(speciesTable_, speciesThermo_),
846 newFormat_(newFormat),
847 imbalanceTol_(ROOTSMALL)
849 read(CHEMKINFileName, thermoFileName, transportFileName);
853 Foam::chemkinReader::chemkinReader
861 speciesTable_(species),
862 reactions_(speciesTable_, speciesThermo_),
868 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
887 chemkinFile = relPath/chemkinFile;
890 if (!thermoFile.empty() && !thermoFile.
isAbsolute())
892 thermoFile = relPath/thermoFile;
897 transportFile = relPath/transportFile;
901 read(chemkinFile, thermoFile, transportFile);