56 const char* Foam::chemkinReader::reactionTypeNames[4] =
60 "nonEquilibriumReversible",
64 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
68 "unimolecularFallOff",
69 "chemicallyActivatedBimolecular",
73 "unknownReactionRateType" 76 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
81 "unknownFallOffFunctionType" 84 void Foam::chemkinReader::initReactionKeywordTable()
86 reactionKeywordTable_.insert(
"M", thirdBodyReactionType);
87 reactionKeywordTable_.insert(
"LOW", unimolecularFallOffReactionType);
88 reactionKeywordTable_.insert
91 chemicallyActivatedBimolecularReactionType
93 reactionKeywordTable_.insert(
"TROE", TroeReactionType);
94 reactionKeywordTable_.insert(
"SRI", SRIReactionType);
95 reactionKeywordTable_.insert(
"LT", LandauTellerReactionType);
96 reactionKeywordTable_.insert(
"RLT", reverseLandauTellerReactionType);
97 reactionKeywordTable_.insert(
"JAN", JanevReactionType);
98 reactionKeywordTable_.insert(
"FIT1", powerSeriesReactionRateType);
99 reactionKeywordTable_.insert(
"HV", radiationActivatedReactionType);
100 reactionKeywordTable_.insert(
"TDEP", speciesTempReactionType);
101 reactionKeywordTable_.insert(
"EXCI", energyLossReactionType);
102 reactionKeywordTable_.insert(
"MOME", plasmaMomentumTransfer);
103 reactionKeywordTable_.insert(
"XSMI", collisionCrossSection);
104 reactionKeywordTable_.insert(
"REV", nonEquilibriumReversibleReactionType);
105 reactionKeywordTable_.insert(
"DUPLICATE", duplicateReactionType);
106 reactionKeywordTable_.insert(
"DUP", duplicateReactionType);
107 reactionKeywordTable_.insert(
"FORD", speciesOrderForward);
108 reactionKeywordTable_.insert(
"RORD", speciesOrderReverse);
109 reactionKeywordTable_.insert(
"UNITS", UnitsOfReaction);
110 reactionKeywordTable_.insert(
"END", end);
114 Foam::scalar Foam::chemkinReader::molecularWeight
116 const List<specieElement>& specieComposition
121 forAll(specieComposition, i)
123 label nAtoms = specieComposition[i].nAtoms();
124 const word& elementName = specieComposition[i].name();
126 if (isotopeAtomicWts_.found(elementName))
128 molWt += nAtoms*isotopeAtomicWts_[elementName];
137 <<
"Unknown element " << elementName
138 <<
" on line " << lineNo_-1 <<
nl 139 <<
" specieComposition: " << specieComposition
148 void Foam::chemkinReader::checkCoeffs
151 const char* reactionRateName,
155 if (reactionCoeffs.size() != nCoeffs)
158 <<
"Wrong number of coefficients for the " << reactionRateName
159 <<
" rate expression on line " 160 << lineNo_-1 <<
", should be " 161 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl 162 <<
"Coefficients are " 163 << reactionCoeffs <<
nl 168 template<
class ReactionRateType>
169 void Foam::chemkinReader::addReactionType
171 const reactionType rType,
172 DynamicList<gasHReaction::specieCoeffs>& lhs,
173 DynamicList<gasHReaction::specieCoeffs>& rhs,
174 const ReactionRateType& rr
183 new IrreversibleReaction
186 ReactionProxy<gasHThermoPhysics>
203 new ReversibleReaction
206 ReactionProxy<gasHThermoPhysics>
224 <<
"Reaction type " << reactionTypeNames[rType]
225 <<
" on line " << lineNo_-1
226 <<
" not handled by this function" 232 <<
"Unknown reaction type " << rType
233 <<
" on line " << lineNo_-1
239 template<
template<
class,
class>
class PressureDependencyType>
240 void Foam::chemkinReader::addPressureDependentReaction
242 const reactionType rType,
243 const fallOffFunctionType fofType,
244 DynamicList<gasHReaction::specieCoeffs>& lhs,
245 DynamicList<gasHReaction::specieCoeffs>& rhs,
249 const HashTable<scalarList>& reactionCoeffsTable,
250 const scalar Afactor0,
251 const scalar AfactorInf,
255 checkCoeffs(k0Coeffs,
"k0", 3);
256 checkCoeffs(kInfCoeffs,
"kInf", 3);
266 PressureDependencyType
267 <ArrheniusReactionRate, LindemannFallOffFunction>
269 ArrheniusReactionRate
271 Afactor0*k0Coeffs[0],
275 ArrheniusReactionRate
277 AfactorInf*kInfCoeffs[0],
281 LindemannFallOffFunction(),
282 thirdBodyEfficiencies(speciesTable_, efficiencies)
291 reactionCoeffsTable[fallOffFunctionNames[fofType]]
294 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
297 <<
"Wrong number of coefficients for Troe rate expression" 298 " on line " << lineNo_-1 <<
", should be 3 or 4 but " 299 << TroeCoeffs.size() <<
" supplied." <<
nl 300 <<
"Coefficients are " 305 if (TroeCoeffs.size() == 3)
307 TroeCoeffs.setSize(4);
308 TroeCoeffs[3] = great;
315 PressureDependencyType
316 <ArrheniusReactionRate, TroeFallOffFunction>
318 ArrheniusReactionRate
320 Afactor0*k0Coeffs[0],
324 ArrheniusReactionRate
326 AfactorInf*kInfCoeffs[0],
337 thirdBodyEfficiencies(speciesTable_, efficiencies)
346 reactionCoeffsTable[fallOffFunctionNames[fofType]]
349 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
352 <<
"Wrong number of coefficients for SRI rate expression" 353 " on line " << lineNo_-1 <<
", should be 3 or 5 but " 354 << SRICoeffs.size() <<
" supplied." <<
nl 355 <<
"Coefficients are " 360 if (SRICoeffs.size() == 3)
362 SRICoeffs.setSize(5);
371 PressureDependencyType
372 <ArrheniusReactionRate, SRIFallOffFunction>
374 ArrheniusReactionRate
376 Afactor0*k0Coeffs[0],
380 ArrheniusReactionRate
382 AfactorInf*kInfCoeffs[0],
394 thirdBodyEfficiencies(speciesTable_, efficiencies)
402 <<
"Fall-off function type " 403 << fallOffFunctionNames[fofType]
404 <<
" on line " << lineNo_-1
405 <<
" not implemented" 412 void Foam::chemkinReader::addReaction
414 DynamicList<gasHReaction::specieCoeffs>& lhs,
415 DynamicList<gasHReaction::specieCoeffs>& rhs,
417 const reactionType rType,
418 const reactionRateType rrType,
419 const fallOffFunctionType fofType,
421 HashTable<scalarList>& reactionCoeffsTable,
425 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
431 const List<specieElement>& specieComposition =
432 speciesComposition_[speciesTable_[lhs[i].index]];
434 forAll(specieComposition, j)
436 label elementi = elementIndices_[specieComposition[j].name()];
438 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
444 const List<specieElement>& specieComposition =
445 speciesComposition_[speciesTable_[rhs[i].index]];
447 forAll(specieComposition, j)
449 label elementi = elementIndices_[specieComposition[j].name()];
451 rhs[i].stoichCoeff*specieComposition[j].nAtoms();
458 const scalar concFactor = 0.001;
462 sumExp += lhs[i].exponent;
464 scalar Afactor =
pow(concFactor, sumExp - 1.0);
466 scalar AfactorRev = Afactor;
468 if (rType == nonEquilibriumReversible)
473 sumExp += rhs[i].exponent;
475 AfactorRev =
pow(concFactor, sumExp - 1.0);
482 if (rType == nonEquilibriumReversible)
485 reactionCoeffsTable[reactionTypeNames[rType]];
487 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
491 new NonEquilibriumReversibleReaction
494 ReactionProxy<gasHThermoPhysics>
501 ArrheniusReactionRate
503 Afactor*ArrheniusCoeffs[0],
505 ArrheniusCoeffs[2]/RR
507 ArrheniusReactionRate
509 AfactorRev*reverseArrheniusCoeffs[0],
510 reverseArrheniusCoeffs[1],
511 reverseArrheniusCoeffs[2]/RR
522 ArrheniusReactionRate
524 Afactor*ArrheniusCoeffs[0],
526 ArrheniusCoeffs[2]/RR
532 case thirdBodyArrhenius:
534 if (rType == nonEquilibriumReversible)
537 reactionCoeffsTable[reactionTypeNames[rType]];
539 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
543 new NonEquilibriumReversibleReaction
547 thirdBodyArrheniusReactionRate
550 ReactionProxy<gasHThermoPhysics>
557 thirdBodyArrheniusReactionRate
559 Afactor*concFactor*ArrheniusCoeffs[0],
561 ArrheniusCoeffs[2]/RR,
562 thirdBodyEfficiencies(speciesTable_, efficiencies)
564 thirdBodyArrheniusReactionRate
566 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
567 reverseArrheniusCoeffs[1],
568 reverseArrheniusCoeffs[2]/RR,
569 thirdBodyEfficiencies(speciesTable_, efficiencies)
580 thirdBodyArrheniusReactionRate
582 Afactor*concFactor*ArrheniusCoeffs[0],
584 ArrheniusCoeffs[2]/RR,
585 thirdBodyEfficiencies(speciesTable_, efficiencies)
591 case unimolecularFallOff:
593 addPressureDependentReaction<FallOffReactionRate>
600 reactionCoeffsTable[reactionRateTypeNames[rrType]],
609 case chemicallyActivatedBimolecular:
611 addPressureDependentReaction<ChemicallyActivatedReactionRate>
619 reactionCoeffsTable[reactionRateTypeNames[rrType]],
630 reactionCoeffsTable[reactionRateTypeNames[rrType]];
631 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
633 if (rType == nonEquilibriumReversible)
636 reactionCoeffsTable[reactionTypeNames[rType]];
637 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
642 word(reactionTypeNames[rType])
643 + reactionRateTypeNames[rrType]
645 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
649 new NonEquilibriumReversibleReaction
653 LandauTellerReactionRate
656 ReactionProxy<gasHThermoPhysics>
663 LandauTellerReactionRate
665 Afactor*ArrheniusCoeffs[0],
667 ArrheniusCoeffs[2]/RR,
668 LandauTellerCoeffs[0],
669 LandauTellerCoeffs[1]
671 LandauTellerReactionRate
673 AfactorRev*reverseArrheniusCoeffs[0],
674 reverseArrheniusCoeffs[1],
675 reverseArrheniusCoeffs[2]/RR,
676 reverseLandauTellerCoeffs[0],
677 reverseLandauTellerCoeffs[1]
688 LandauTellerReactionRate
690 Afactor*ArrheniusCoeffs[0],
692 ArrheniusCoeffs[2]/RR,
693 LandauTellerCoeffs[0],
694 LandauTellerCoeffs[1]
703 reactionCoeffsTable[reactionRateTypeNames[rrType]];
705 checkCoeffs(JanevCoeffs,
"Janev", 9);
713 Afactor*ArrheniusCoeffs[0],
715 ArrheniusCoeffs[2]/RR,
716 FixedList<scalar, 9>(JanevCoeffs)
724 reactionCoeffsTable[reactionRateTypeNames[rrType]];
726 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
732 powerSeriesReactionRate
734 Afactor*ArrheniusCoeffs[0],
736 ArrheniusCoeffs[2]/RR,
737 FixedList<scalar, 4>(powerSeriesCoeffs)
742 case unknownReactionRateType:
745 <<
"Internal error on line " << lineNo_-1
746 <<
": reaction rate type has not been set" 753 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
754 <<
" on line " << lineNo_-1
755 <<
" not implemented" 763 if (
mag(nAtoms[i]) > imbalanceTol_)
766 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
767 <<
" in " << elementNames_[i]
768 <<
" in reaction" <<
nl 769 << reactions_.last() <<
nl 770 <<
" on line " << lineNo_-1
777 reactionCoeffsTable.clear();
781 void Foam::chemkinReader::read
783 const fileName& CHEMKINFileName,
784 const fileName& thermoFileName,
785 const fileName& transportFileName
791 transportDict_.read(IFstream(transportFileName)());
795 std::ifstream thermoStream(thermoFileName.c_str());
800 <<
"file " << thermoFileName <<
" not found" 804 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
805 yy_switch_to_buffer(bufferPtr);
810 yy_delete_buffer(bufferPtr);
815 std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
820 <<
"file " << CHEMKINFileName <<
" not found" 824 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
825 yy_switch_to_buffer(bufferPtr);
827 initReactionKeywordTable();
832 yy_delete_buffer(bufferPtr);
838 Foam::chemkinReader::chemkinReader
849 speciesTable_(species),
850 reactions_(speciesTable_, speciesThermo_),
851 newFormat_(newFormat),
852 imbalanceTol_(rootSmall)
854 read(CHEMKINFileName, thermoFileName, transportFileName);
858 Foam::chemkinReader::chemkinReader
866 speciesTable_(species),
867 reactions_(speciesTable_, speciesThermo_),
869 imbalanceTol_(thermoDict.
lookupOrDefault(
"imbalanceTolerance", rootSmall))
873 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
880 if (thermoDict.
found(
"CHEMKINThermoFile"))
894 if (!chemkinFile.isAbsolute())
896 chemkinFile = relPath/chemkinFile;
901 thermoFile = relPath/thermoFile;
904 if (!transportFile.isAbsolute())
906 transportFile = relPath/transportFile;
910 read(chemkinFile, thermoFile, transportFile);
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
#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.
A class for handling file names.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A list of keyword definitions, which are a keyword followed by any number of values (e...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static const fileName null
An empty fileName.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Macros for easy insertion into run-time selection tables.
static scalar TlowDefault
Default temperature limits of applicability of reaction rates.
bool read(const char *, int32_t &)
const fileName & name() const
Return the dictionary name.
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.
addChemistryReaderType(chemkinReader, gasHThermoPhysics)
static scalar ThighDefault
atomicWeightTable atomicWeights
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
A wordList with hashed indices for faster lookup by name.
dimensioned< scalar > mag(const dimensioned< Type > &)
fileName path() const
Return directory path name (part before last /)
const scalar RR
Universal gas constant (default in [J/(kmol K)])
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
sutherlandTransport< species::thermo< janafThermo< perfectGas< specie > >, sensibleEnthalpy > > gasHThermoPhysics