57 const char* Foam::chemkinReader::reactionTypeNames[4] =
61 "nonEquilibriumReversible",
65 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
69 "unimolecularFallOff",
70 "chemicallyActivatedBimolecular",
74 "unknownReactionRateType"
77 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
82 "unknownFallOffFunctionType"
85 void Foam::chemkinReader::initReactionKeywordTable()
87 reactionKeywordTable_.
insert(
"M", thirdBodyReactionType);
88 reactionKeywordTable_.
insert(
"LOW", unimolecularFallOffReactionType);
89 reactionKeywordTable_.
insert
92 chemicallyActivatedBimolecularReactionType
94 reactionKeywordTable_.
insert(
"TROE", TroeReactionType);
95 reactionKeywordTable_.
insert(
"SRI", SRIReactionType);
96 reactionKeywordTable_.
insert(
"LT", LandauTellerReactionType);
97 reactionKeywordTable_.
insert(
"RLT", reverseLandauTellerReactionType);
98 reactionKeywordTable_.
insert(
"JAN", JanevReactionType);
99 reactionKeywordTable_.
insert(
"FIT1", powerSeriesReactionRateType);
100 reactionKeywordTable_.
insert(
"HV", radiationActivatedReactionType);
101 reactionKeywordTable_.
insert(
"TDEP", speciesTempReactionType);
102 reactionKeywordTable_.
insert(
"EXCI", energyLossReactionType);
103 reactionKeywordTable_.
insert(
"MOME", plasmaMomentumTransfer);
104 reactionKeywordTable_.
insert(
"XSMI", collisionCrossSection);
105 reactionKeywordTable_.
insert(
"REV", nonEquilibriumReversibleReactionType);
106 reactionKeywordTable_.
insert(
"DUPLICATE", duplicateReactionType);
107 reactionKeywordTable_.
insert(
"DUP", duplicateReactionType);
108 reactionKeywordTable_.
insert(
"FORD", speciesOrderForward);
109 reactionKeywordTable_.
insert(
"RORD", speciesOrderReverse);
110 reactionKeywordTable_.
insert(
"UNITS", UnitsOfReaction);
111 reactionKeywordTable_.
insert(
"END", end);
115 Foam::scalar Foam::chemkinReader::molecularWeight
117 const List<specieElement>& specieComposition
122 forAll(specieComposition, i)
124 label nAtoms = specieComposition[i].nAtoms();
125 const word& elementName = specieComposition[i].name();
127 if (isotopeAtomicWts_.found(elementName))
129 molWt += nAtoms*isotopeAtomicWts_[elementName];
138 <<
"Unknown element " << elementName
139 <<
" on line " << lineNo_-1 <<
nl
140 <<
" specieComposition: " << specieComposition
149 void Foam::chemkinReader::checkCoeffs
152 const char* reactionRateName,
156 if (reactionCoeffs.size() != nCoeffs)
159 <<
"Wrong number of coefficients for the " << reactionRateName
160 <<
" rate expression on line "
161 << lineNo_-1 <<
", should be "
162 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl
163 <<
"Coefficients are "
164 << reactionCoeffs <<
nl
169 template<
class ReactionRateType>
170 void Foam::chemkinReader::addReactionType
172 const reactionType rType,
173 DynamicList<gasHReaction::specieCoeffs>& lhs,
174 DynamicList<gasHReaction::specieCoeffs>& rhs,
175 const ReactionRateType& rr
184 new IrreversibleReaction
187 Reaction<gasHThermoPhysics>
204 new ReversibleReaction
207 Reaction<gasHThermoPhysics>
225 <<
"Reaction type " << reactionTypeNames[rType]
226 <<
" on line " << lineNo_-1
227 <<
" not handled by this function"
233 <<
"Unknown reaction type " << rType
234 <<
" on line " << lineNo_-1
240 template<
template<
class,
class>
class PressureDependencyType>
241 void Foam::chemkinReader::addPressureDependentReaction
243 const reactionType rType,
244 const fallOffFunctionType fofType,
245 DynamicList<gasHReaction::specieCoeffs>& lhs,
246 DynamicList<gasHReaction::specieCoeffs>& rhs,
250 const HashTable<scalarList>& reactionCoeffsTable,
251 const scalar Afactor0,
252 const scalar AfactorInf,
256 checkCoeffs(k0Coeffs,
"k0", 3);
257 checkCoeffs(kInfCoeffs,
"kInf", 3);
267 PressureDependencyType
268 <ArrheniusReactionRate, LindemannFallOffFunction>
270 ArrheniusReactionRate
272 Afactor0*k0Coeffs[0],
276 ArrheniusReactionRate
278 AfactorInf*kInfCoeffs[0],
282 LindemannFallOffFunction(),
283 thirdBodyEfficiencies(speciesTable_, efficiencies)
292 reactionCoeffsTable[fallOffFunctionNames[fofType]]
295 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
298 <<
"Wrong number of coefficients for Troe rate expression"
299 " on line " << lineNo_-1 <<
", should be 3 or 4 but "
300 << TroeCoeffs.size() <<
" supplied." <<
nl
301 <<
"Coefficients are "
306 if (TroeCoeffs.size() == 3)
308 TroeCoeffs.setSize(4);
309 TroeCoeffs[3] = GREAT;
316 PressureDependencyType
317 <ArrheniusReactionRate, TroeFallOffFunction>
319 ArrheniusReactionRate
321 Afactor0*k0Coeffs[0],
325 ArrheniusReactionRate
327 AfactorInf*kInfCoeffs[0],
338 thirdBodyEfficiencies(speciesTable_, efficiencies)
347 reactionCoeffsTable[fallOffFunctionNames[fofType]]
350 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
353 <<
"Wrong number of coefficients for SRI rate expression"
354 " on line " << lineNo_-1 <<
", should be 3 or 5 but "
355 << SRICoeffs.size() <<
" supplied." <<
nl
356 <<
"Coefficients are "
361 if (SRICoeffs.size() == 3)
363 SRICoeffs.setSize(5);
372 PressureDependencyType
373 <ArrheniusReactionRate, SRIFallOffFunction>
375 ArrheniusReactionRate
377 Afactor0*k0Coeffs[0],
381 ArrheniusReactionRate
383 AfactorInf*kInfCoeffs[0],
395 thirdBodyEfficiencies(speciesTable_, efficiencies)
403 <<
"Fall-off function type "
404 << fallOffFunctionNames[fofType]
405 <<
" on line " << lineNo_-1
406 <<
" not implemented"
413 void Foam::chemkinReader::addReaction
415 DynamicList<gasHReaction::specieCoeffs>& lhs,
416 DynamicList<gasHReaction::specieCoeffs>& rhs,
418 const reactionType rType,
419 const reactionRateType rrType,
420 const fallOffFunctionType fofType,
422 HashTable<scalarList>& reactionCoeffsTable,
426 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
432 const List<specieElement>& specieComposition =
433 speciesComposition_[speciesTable_[lhs[i].index]];
435 forAll(specieComposition, j)
437 label elementi = elementIndices_[specieComposition[j].name()];
439 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
445 const List<specieElement>& specieComposition =
446 speciesComposition_[speciesTable_[rhs[i].index]];
448 forAll(specieComposition, j)
450 label elementi = elementIndices_[specieComposition[j].name()];
452 rhs[i].stoichCoeff*specieComposition[j].nAtoms();
459 const scalar concFactor = 0.001;
463 sumExp += lhs[i].exponent;
465 scalar Afactor =
pow(concFactor, sumExp - 1.0);
467 scalar AfactorRev = Afactor;
469 if (rType == nonEquilibriumReversible)
474 sumExp += rhs[i].exponent;
476 AfactorRev =
pow(concFactor, sumExp - 1.0);
483 if (rType == nonEquilibriumReversible)
486 reactionCoeffsTable[reactionTypeNames[rType]];
488 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
492 new NonEquilibriumReversibleReaction
495 Reaction<gasHThermoPhysics>
502 ArrheniusReactionRate
504 Afactor*ArrheniusCoeffs[0],
506 ArrheniusCoeffs[2]/
RR
508 ArrheniusReactionRate
510 AfactorRev*reverseArrheniusCoeffs[0],
511 reverseArrheniusCoeffs[1],
512 reverseArrheniusCoeffs[2]/
RR
523 ArrheniusReactionRate
525 Afactor*ArrheniusCoeffs[0],
527 ArrheniusCoeffs[2]/
RR
533 case thirdBodyArrhenius:
535 if (rType == nonEquilibriumReversible)
538 reactionCoeffsTable[reactionTypeNames[rType]];
540 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
544 new NonEquilibriumReversibleReaction
548 thirdBodyArrheniusReactionRate
551 Reaction<gasHThermoPhysics>
558 thirdBodyArrheniusReactionRate
560 Afactor*concFactor*ArrheniusCoeffs[0],
562 ArrheniusCoeffs[2]/
RR,
563 thirdBodyEfficiencies(speciesTable_, efficiencies)
565 thirdBodyArrheniusReactionRate
567 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
568 reverseArrheniusCoeffs[1],
569 reverseArrheniusCoeffs[2]/
RR,
570 thirdBodyEfficiencies(speciesTable_, efficiencies)
581 thirdBodyArrheniusReactionRate
583 Afactor*concFactor*ArrheniusCoeffs[0],
585 ArrheniusCoeffs[2]/
RR,
586 thirdBodyEfficiencies(speciesTable_, efficiencies)
592 case unimolecularFallOff:
594 addPressureDependentReaction<FallOffReactionRate>
601 reactionCoeffsTable[reactionRateTypeNames[rrType]],
610 case chemicallyActivatedBimolecular:
612 addPressureDependentReaction<ChemicallyActivatedReactionRate>
620 reactionCoeffsTable[reactionRateTypeNames[rrType]],
631 reactionCoeffsTable[reactionRateTypeNames[rrType]];
632 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
634 if (rType == nonEquilibriumReversible)
637 reactionCoeffsTable[reactionTypeNames[rType]];
638 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
643 word(reactionTypeNames[rType])
644 + reactionRateTypeNames[rrType]
646 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
650 new NonEquilibriumReversibleReaction
653 Reaction<gasHThermoPhysics>
660 LandauTellerReactionRate
662 Afactor*ArrheniusCoeffs[0],
664 ArrheniusCoeffs[2]/
RR,
665 LandauTellerCoeffs[0],
666 LandauTellerCoeffs[1]
668 LandauTellerReactionRate
670 AfactorRev*reverseArrheniusCoeffs[0],
671 reverseArrheniusCoeffs[1],
672 reverseArrheniusCoeffs[2]/
RR,
673 reverseLandauTellerCoeffs[0],
674 reverseLandauTellerCoeffs[1]
685 LandauTellerReactionRate
687 Afactor*ArrheniusCoeffs[0],
689 ArrheniusCoeffs[2]/
RR,
690 LandauTellerCoeffs[0],
691 LandauTellerCoeffs[1]
700 reactionCoeffsTable[reactionRateTypeNames[rrType]];
702 checkCoeffs(JanevCoeffs,
"Janev", 9);
710 Afactor*ArrheniusCoeffs[0],
712 ArrheniusCoeffs[2]/
RR,
713 FixedList<scalar, 9>(JanevCoeffs)
721 reactionCoeffsTable[reactionRateTypeNames[rrType]];
723 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
729 powerSeriesReactionRate
731 Afactor*ArrheniusCoeffs[0],
733 ArrheniusCoeffs[2]/
RR,
734 FixedList<scalar, 4>(powerSeriesCoeffs)
739 case unknownReactionRateType:
742 <<
"Internal error on line " << lineNo_-1
743 <<
": reaction rate type has not been set"
750 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
751 <<
" on line " << lineNo_-1
752 <<
" not implemented"
760 if (
mag(nAtoms[i]) > imbalanceTol_)
763 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
764 <<
" in " << elementNames_[i]
765 <<
" in reaction" <<
nl
766 << reactions_.last() <<
nl
767 <<
" on line " << lineNo_-1
774 reactionCoeffsTable.clear();
778 void Foam::chemkinReader::read
780 const fileName& CHEMKINFileName,
781 const fileName& thermoFileName,
782 const fileName& transportFileName
785 transportDict_.read(IFstream(transportFileName)());
787 if (!thermoFileName.empty())
789 std::ifstream thermoStream(thermoFileName);
794 <<
"file " << thermoFileName <<
" not found"
798 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
799 yy_switch_to_buffer(bufferPtr);
804 yy_delete_buffer(bufferPtr);
809 std::ifstream CHEMKINStream(CHEMKINFileName);
814 <<
"file " << CHEMKINFileName <<
" not found"
818 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
819 yy_switch_to_buffer(bufferPtr);
821 initReactionKeywordTable();
826 yy_delete_buffer(bufferPtr);
832 Foam::chemkinReader::chemkinReader
843 speciesTable_(species),
844 reactions_(speciesTable_, speciesThermo_),
845 newFormat_(newFormat),
846 imbalanceTol_(ROOTSMALL)
848 read(CHEMKINFileName, thermoFileName, transportFileName);
852 Foam::chemkinReader::chemkinReader
860 speciesTable_(species),
861 reactions_(speciesTable_, speciesThermo_),
867 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
886 chemkinFile = relPath/chemkinFile;
891 thermoFile = relPath/thermoFile;
896 transportFile = relPath/transportFile;
900 read(chemkinFile, thermoFile, transportFile);