55 const char* Foam::chemkinReader::reactionTypeNames[4] =
59 "nonEquilibriumReversible",
63 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
67 "unimolecularFallOff",
68 "chemicallyActivatedBimolecular",
72 "unknownReactionRateType" 75 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
80 "unknownFallOffFunctionType" 83 void Foam::chemkinReader::initReactionKeywordTable()
85 reactionKeywordTable_.insert(
"M", thirdBodyReactionType);
86 reactionKeywordTable_.insert(
"LOW", unimolecularFallOffReactionType);
87 reactionKeywordTable_.insert
90 chemicallyActivatedBimolecularReactionType
92 reactionKeywordTable_.insert(
"TROE", TroeReactionType);
93 reactionKeywordTable_.insert(
"SRI", SRIReactionType);
94 reactionKeywordTable_.insert(
"LT", LandauTellerReactionType);
95 reactionKeywordTable_.insert(
"RLT", reverseLandauTellerReactionType);
96 reactionKeywordTable_.insert(
"JAN", JanevReactionType);
97 reactionKeywordTable_.insert(
"FIT1", powerSeriesReactionRateType);
98 reactionKeywordTable_.insert(
"HV", radiationActivatedReactionType);
99 reactionKeywordTable_.insert(
"TDEP", speciesTempReactionType);
100 reactionKeywordTable_.insert(
"EXCI", energyLossReactionType);
101 reactionKeywordTable_.insert(
"MOME", plasmaMomentumTransfer);
102 reactionKeywordTable_.insert(
"XSMI", collisionCrossSection);
103 reactionKeywordTable_.insert(
"REV", nonEquilibriumReversibleReactionType);
104 reactionKeywordTable_.insert(
"DUPLICATE", duplicateReactionType);
105 reactionKeywordTable_.insert(
"DUP", duplicateReactionType);
106 reactionKeywordTable_.insert(
"FORD", speciesOrderForward);
107 reactionKeywordTable_.insert(
"RORD", speciesOrderReverse);
108 reactionKeywordTable_.insert(
"UNITS", UnitsOfReaction);
109 reactionKeywordTable_.insert(
"END", end);
113 Foam::scalar Foam::chemkinReader::molecularWeight
115 const List<specieElement>& specieComposition
120 forAll(specieComposition, i)
122 label nAtoms = specieComposition[i].nAtoms();
123 const word& elementName = specieComposition[i].name();
125 if (isotopeAtomicWts_.found(elementName))
127 molWt += nAtoms*isotopeAtomicWts_[elementName];
136 <<
"Unknown element " << elementName
137 <<
" on line " << lineNo_-1 <<
nl 138 <<
" specieComposition: " << specieComposition
147 void Foam::chemkinReader::checkCoeffs
150 const char* reactionRateName,
154 if (reactionCoeffs.size() != nCoeffs)
157 <<
"Wrong number of coefficients for the " << reactionRateName
158 <<
" rate expression on line " 159 << lineNo_-1 <<
", should be " 160 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl 161 <<
"Coefficients are " 162 << reactionCoeffs <<
nl 167 template<
class ReactionRateType>
168 void Foam::chemkinReader::addReactionType
170 const reactionType rType,
171 DynamicList<gasHReaction::specieCoeffs>& lhs,
172 DynamicList<gasHReaction::specieCoeffs>& rhs,
173 const ReactionRateType& rr
182 new IrreversibleReaction
185 Reaction<gasHThermoPhysics>
202 new ReversibleReaction
205 Reaction<gasHThermoPhysics>
223 <<
"Reaction type " << reactionTypeNames[rType]
224 <<
" on line " << lineNo_-1
225 <<
" not handled by this function" 231 <<
"Unknown reaction type " << rType
232 <<
" on line " << lineNo_-1
238 template<
template<
class,
class>
class PressureDependencyType>
239 void Foam::chemkinReader::addPressureDependentReaction
241 const reactionType rType,
242 const fallOffFunctionType fofType,
243 DynamicList<gasHReaction::specieCoeffs>& lhs,
244 DynamicList<gasHReaction::specieCoeffs>& rhs,
248 const HashTable<scalarList>& reactionCoeffsTable,
249 const scalar Afactor0,
250 const scalar AfactorInf,
254 checkCoeffs(k0Coeffs,
"k0", 3);
255 checkCoeffs(kInfCoeffs,
"kInf", 3);
265 PressureDependencyType
266 <ArrheniusReactionRate, LindemannFallOffFunction>
268 ArrheniusReactionRate
270 Afactor0*k0Coeffs[0],
274 ArrheniusReactionRate
276 AfactorInf*kInfCoeffs[0],
280 LindemannFallOffFunction(),
281 thirdBodyEfficiencies(speciesTable_, efficiencies)
290 reactionCoeffsTable[fallOffFunctionNames[fofType]]
293 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
296 <<
"Wrong number of coefficients for Troe rate expression" 297 " on line " << lineNo_-1 <<
", should be 3 or 4 but " 298 << TroeCoeffs.size() <<
" supplied." <<
nl 299 <<
"Coefficients are " 304 if (TroeCoeffs.size() == 3)
306 TroeCoeffs.setSize(4);
307 TroeCoeffs[3] = GREAT;
314 PressureDependencyType
315 <ArrheniusReactionRate, TroeFallOffFunction>
317 ArrheniusReactionRate
319 Afactor0*k0Coeffs[0],
323 ArrheniusReactionRate
325 AfactorInf*kInfCoeffs[0],
336 thirdBodyEfficiencies(speciesTable_, efficiencies)
345 reactionCoeffsTable[fallOffFunctionNames[fofType]]
348 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
351 <<
"Wrong number of coefficients for SRI rate expression" 352 " on line " << lineNo_-1 <<
", should be 3 or 5 but " 353 << SRICoeffs.size() <<
" supplied." <<
nl 354 <<
"Coefficients are " 359 if (SRICoeffs.size() == 3)
361 SRICoeffs.setSize(5);
370 PressureDependencyType
371 <ArrheniusReactionRate, SRIFallOffFunction>
373 ArrheniusReactionRate
375 Afactor0*k0Coeffs[0],
379 ArrheniusReactionRate
381 AfactorInf*kInfCoeffs[0],
393 thirdBodyEfficiencies(speciesTable_, efficiencies)
401 <<
"Fall-off function type " 402 << fallOffFunctionNames[fofType]
403 <<
" on line " << lineNo_-1
404 <<
" not implemented" 411 void Foam::chemkinReader::addReaction
413 DynamicList<gasHReaction::specieCoeffs>& lhs,
414 DynamicList<gasHReaction::specieCoeffs>& rhs,
416 const reactionType rType,
417 const reactionRateType rrType,
418 const fallOffFunctionType fofType,
420 HashTable<scalarList>& reactionCoeffsTable,
424 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
430 const List<specieElement>& specieComposition =
431 speciesComposition_[speciesTable_[lhs[i].index]];
433 forAll(specieComposition, j)
435 label elementi = elementIndices_[specieComposition[j].name()];
437 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
443 const List<specieElement>& specieComposition =
444 speciesComposition_[speciesTable_[rhs[i].index]];
446 forAll(specieComposition, j)
448 label elementi = elementIndices_[specieComposition[j].name()];
450 rhs[i].stoichCoeff*specieComposition[j].nAtoms();
457 const scalar concFactor = 0.001;
461 sumExp += lhs[i].exponent;
463 scalar Afactor =
pow(concFactor, sumExp - 1.0);
465 scalar AfactorRev = Afactor;
467 if (rType == nonEquilibriumReversible)
472 sumExp += rhs[i].exponent;
474 AfactorRev =
pow(concFactor, sumExp - 1.0);
481 if (rType == nonEquilibriumReversible)
484 reactionCoeffsTable[reactionTypeNames[rType]];
486 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
490 new NonEquilibriumReversibleReaction
493 Reaction<gasHThermoPhysics>
500 ArrheniusReactionRate
502 Afactor*ArrheniusCoeffs[0],
504 ArrheniusCoeffs[2]/RR
506 ArrheniusReactionRate
508 AfactorRev*reverseArrheniusCoeffs[0],
509 reverseArrheniusCoeffs[1],
510 reverseArrheniusCoeffs[2]/RR
521 ArrheniusReactionRate
523 Afactor*ArrheniusCoeffs[0],
525 ArrheniusCoeffs[2]/RR
531 case thirdBodyArrhenius:
533 if (rType == nonEquilibriumReversible)
536 reactionCoeffsTable[reactionTypeNames[rType]];
538 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
542 new NonEquilibriumReversibleReaction
546 thirdBodyArrheniusReactionRate
549 Reaction<gasHThermoPhysics>
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)
590 case unimolecularFallOff:
592 addPressureDependentReaction<FallOffReactionRate>
599 reactionCoeffsTable[reactionRateTypeNames[rrType]],
608 case chemicallyActivatedBimolecular:
610 addPressureDependentReaction<ChemicallyActivatedReactionRate>
618 reactionCoeffsTable[reactionRateTypeNames[rrType]],
629 reactionCoeffsTable[reactionRateTypeNames[rrType]];
630 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
632 if (rType == nonEquilibriumReversible)
635 reactionCoeffsTable[reactionTypeNames[rType]];
636 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
641 word(reactionTypeNames[rType])
642 + reactionRateTypeNames[rrType]
644 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
648 new NonEquilibriumReversibleReaction
651 Reaction<gasHThermoPhysics>
658 LandauTellerReactionRate
660 Afactor*ArrheniusCoeffs[0],
662 ArrheniusCoeffs[2]/RR,
663 LandauTellerCoeffs[0],
664 LandauTellerCoeffs[1]
666 LandauTellerReactionRate
668 AfactorRev*reverseArrheniusCoeffs[0],
669 reverseArrheniusCoeffs[1],
670 reverseArrheniusCoeffs[2]/RR,
671 reverseLandauTellerCoeffs[0],
672 reverseLandauTellerCoeffs[1]
683 LandauTellerReactionRate
685 Afactor*ArrheniusCoeffs[0],
687 ArrheniusCoeffs[2]/RR,
688 LandauTellerCoeffs[0],
689 LandauTellerCoeffs[1]
698 reactionCoeffsTable[reactionRateTypeNames[rrType]];
700 checkCoeffs(JanevCoeffs,
"Janev", 9);
708 Afactor*ArrheniusCoeffs[0],
710 ArrheniusCoeffs[2]/RR,
711 FixedList<scalar, 9>(JanevCoeffs)
719 reactionCoeffsTable[reactionRateTypeNames[rrType]];
721 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
727 powerSeriesReactionRate
729 Afactor*ArrheniusCoeffs[0],
731 ArrheniusCoeffs[2]/RR,
732 FixedList<scalar, 4>(powerSeriesCoeffs)
737 case unknownReactionRateType:
740 <<
"Internal error on line " << lineNo_-1
741 <<
": reaction rate type has not been set" 748 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
749 <<
" on line " << lineNo_-1
750 <<
" not implemented" 758 if (
mag(nAtoms[i]) > imbalanceTol_)
761 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
762 <<
" in " << elementNames_[i]
763 <<
" in reaction" <<
nl 764 << reactions_.last() <<
nl 765 <<
" on line " << lineNo_-1
772 reactionCoeffsTable.clear();
776 void Foam::chemkinReader::read
778 const fileName& CHEMKINFileName,
779 const fileName& thermoFileName,
780 const fileName& transportFileName
783 transportDict_.read(IFstream(transportFileName)());
787 std::ifstream thermoStream(thermoFileName.c_str());
792 <<
"file " << thermoFileName <<
" not found" 796 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
797 yy_switch_to_buffer(bufferPtr);
802 yy_delete_buffer(bufferPtr);
807 std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
812 <<
"file " << CHEMKINFileName <<
" not found" 816 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
817 yy_switch_to_buffer(bufferPtr);
819 initReactionKeywordTable();
824 yy_delete_buffer(bufferPtr);
830 Foam::chemkinReader::chemkinReader
841 speciesTable_(species),
842 reactions_(speciesTable_, speciesThermo_),
843 newFormat_(newFormat),
844 imbalanceTol_(ROOTSMALL)
846 read(CHEMKINFileName, thermoFileName, transportFileName);
850 Foam::chemkinReader::chemkinReader
858 speciesTable_(species),
859 reactions_(speciesTable_, speciesThermo_),
861 imbalanceTol_(thermoDict.
lookupOrDefault(
"imbalanceTolerance", ROOTSMALL))
865 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
872 if (thermoDict.
found(
"CHEMKINThermoFile"))
886 if (!chemkinFile.isAbsolute())
888 chemkinFile = relPath/chemkinFile;
893 thermoFile = relPath/thermoFile;
896 if (!transportFile.isAbsolute())
898 transportFile = relPath/transportFile;
902 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 occurences 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.
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)
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