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 ArrheniusCoeffs[2]/RR,
728 FixedList<scalar, 4>(powerSeriesCoeffs)
733 case unknownReactionRateType:
736 <<
"Internal error on line " << lineNo_-1
737 <<
": reaction rate type has not been set" 744 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
745 <<
" on line " << lineNo_-1
746 <<
" not implemented" 754 if (
mag(nAtoms[i]) > imbalanceTol_)
757 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
758 <<
" in " << elementNames_[i]
759 <<
" in reaction" <<
nl 760 << reactions_.last() <<
nl 761 <<
" on line " << lineNo_-1
768 reactionCoeffsTable.clear();
772 void Foam::chemkinReader::read
774 const fileName& CHEMKINFileName,
775 const fileName& thermoFileName,
776 const fileName& transportFileName
782 transportDict_.
read(IFstream(transportFileName)());
786 std::ifstream thermoStream(thermoFileName.c_str());
791 <<
"file " << thermoFileName <<
" not found" 795 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
796 yy_switch_to_buffer(bufferPtr);
801 yy_delete_buffer(bufferPtr);
806 std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
811 <<
"file " << CHEMKINFileName <<
" not found" 815 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
816 yy_switch_to_buffer(bufferPtr);
818 initReactionKeywordTable();
823 yy_delete_buffer(bufferPtr);
829 Foam::chemkinReader::chemkinReader
832 const fileName& CHEMKINFileName,
833 const fileName& transportFileName,
834 const fileName& thermoFileName,
840 speciesTable_(species),
841 newFormat_(newFormat),
842 imbalanceTol_(rootSmall)
844 read(CHEMKINFileName, thermoFileName, transportFileName);
848 Foam::chemkinReader::chemkinReader
850 const dictionary& thermoDict,
856 speciesTable_(species),
857 newFormat_(thermoDict.lookupOrDefault(
"newFormat", false)),
858 imbalanceTol_(thermoDict.lookupOrDefault(
"imbalanceTolerance", rootSmall))
862 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
865 fileName chemkinFile(fileName(thermoDict.lookup(
"CHEMKINFile")).
expand());
869 if (thermoDict.found(
"CHEMKINThermoFile"))
871 thermoFile = fileName(thermoDict.lookup(
"CHEMKINThermoFile")).
expand();
874 fileName transportFile
876 fileName(thermoDict.lookup(
"CHEMKINTransportFile")).
expand()
880 fileName relPath = thermoDict.name().path();
883 if (!chemkinFile.isAbsolute())
885 chemkinFile = relPath/chemkinFile;
890 thermoFile = relPath/thermoFile;
893 if (!transportFile.isAbsolute())
895 transportFile = relPath/transportFile;
899 read(chemkinFile, thermoFile, transportFile);
#define forAll(list, i)
Loop across all elements in list.
string expand(const string &, const HashTable< string, word, string::hash > &mapping, const char sigil='$')
Expand occurrences of variables according to the mapping.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void size(const label)
Override size to be inconsistent with allocated storage.
static const fileName null
An empty fileName.
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool read(Istream &, const bool keepHeader=false)
Read dictionary from Istream, optionally keeping the header.
hashedWordList speciesTable
A table of species as a hashedWordList.
Macros for easy insertion into run-time selection tables.
static scalar TlowDefault
Default temperature limits of applicability of reaction rates.
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
bool read(const char *, int32_t &)
bool found(const Key &) const
Return true if hashedEntry is found in table.
bool isAbsolute() const
Return true if file name is absolute.
List< scalar > scalarList
A List of scalars.
static scalar ThighDefault
sutherlandTransport< species::thermo< janafThermo< perfectGas< specie > >, sensibleEnthalpy > > thermoPhysics
atomicWeightTable atomicWeights
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
const scalar RR
Universal gas constant (default in [J/kmol/K])