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].elementName;
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 specieComposition_[speciesTable_[lhs[i].index]];
433 forAll(specieComposition, j)
435 label elementi = elementIndices_[specieComposition[j].elementName];
436 nAtoms[elementi] += lhs[i].stoichCoeff*specieComposition[j].nAtoms;
442 const List<specieElement>& specieComposition =
443 specieComposition_[speciesTable_[rhs[i].index]];
445 forAll(specieComposition, j)
447 label elementi = elementIndices_[specieComposition[j].elementName];
448 nAtoms[elementi] -= rhs[i].stoichCoeff*specieComposition[j].nAtoms;
455 const scalar concFactor = 0.001;
459 sumExp += lhs[i].exponent;
461 scalar Afactor =
pow(concFactor, sumExp - 1.0);
463 scalar AfactorRev = Afactor;
465 if (rType == nonEquilibriumReversible)
470 sumExp += rhs[i].exponent;
472 AfactorRev =
pow(concFactor, sumExp - 1.0);
479 if (rType == nonEquilibriumReversible)
482 reactionCoeffsTable[reactionTypeNames[rType]];
484 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
488 new NonEquilibriumReversibleReaction
491 Reaction<gasHThermoPhysics>
498 ArrheniusReactionRate
500 Afactor*ArrheniusCoeffs[0],
502 ArrheniusCoeffs[2]/RR
504 ArrheniusReactionRate
506 AfactorRev*reverseArrheniusCoeffs[0],
507 reverseArrheniusCoeffs[1],
508 reverseArrheniusCoeffs[2]/RR
519 ArrheniusReactionRate
521 Afactor*ArrheniusCoeffs[0],
523 ArrheniusCoeffs[2]/RR
529 case thirdBodyArrhenius:
531 if (rType == nonEquilibriumReversible)
534 reactionCoeffsTable[reactionTypeNames[rType]];
536 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
540 new NonEquilibriumReversibleReaction
544 thirdBodyArrheniusReactionRate
547 Reaction<gasHThermoPhysics>
554 thirdBodyArrheniusReactionRate
556 Afactor*concFactor*ArrheniusCoeffs[0],
558 ArrheniusCoeffs[2]/RR,
559 thirdBodyEfficiencies(speciesTable_, efficiencies)
561 thirdBodyArrheniusReactionRate
563 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
564 reverseArrheniusCoeffs[1],
565 reverseArrheniusCoeffs[2]/RR,
566 thirdBodyEfficiencies(speciesTable_, efficiencies)
577 thirdBodyArrheniusReactionRate
579 Afactor*concFactor*ArrheniusCoeffs[0],
581 ArrheniusCoeffs[2]/RR,
582 thirdBodyEfficiencies(speciesTable_, efficiencies)
588 case unimolecularFallOff:
590 addPressureDependentReaction<FallOffReactionRate>
597 reactionCoeffsTable[reactionRateTypeNames[rrType]],
606 case chemicallyActivatedBimolecular:
608 addPressureDependentReaction<ChemicallyActivatedReactionRate>
616 reactionCoeffsTable[reactionRateTypeNames[rrType]],
627 reactionCoeffsTable[reactionRateTypeNames[rrType]];
628 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
630 if (rType == nonEquilibriumReversible)
633 reactionCoeffsTable[reactionTypeNames[rType]];
634 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
639 word(reactionTypeNames[rType])
640 + reactionRateTypeNames[rrType]
642 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
646 new NonEquilibriumReversibleReaction
649 Reaction<gasHThermoPhysics>
656 LandauTellerReactionRate
658 Afactor*ArrheniusCoeffs[0],
660 ArrheniusCoeffs[2]/RR,
661 LandauTellerCoeffs[0],
662 LandauTellerCoeffs[1]
664 LandauTellerReactionRate
666 AfactorRev*reverseArrheniusCoeffs[0],
667 reverseArrheniusCoeffs[1],
668 reverseArrheniusCoeffs[2]/RR,
669 reverseLandauTellerCoeffs[0],
670 reverseLandauTellerCoeffs[1]
681 LandauTellerReactionRate
683 Afactor*ArrheniusCoeffs[0],
685 ArrheniusCoeffs[2]/RR,
686 LandauTellerCoeffs[0],
687 LandauTellerCoeffs[1]
696 reactionCoeffsTable[reactionRateTypeNames[rrType]];
698 checkCoeffs(JanevCoeffs,
"Janev", 9);
706 Afactor*ArrheniusCoeffs[0],
708 ArrheniusCoeffs[2]/RR,
709 FixedList<scalar, 9>(JanevCoeffs)
717 reactionCoeffsTable[reactionRateTypeNames[rrType]];
719 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
725 powerSeriesReactionRate
727 Afactor*ArrheniusCoeffs[0],
729 ArrheniusCoeffs[2]/RR,
730 FixedList<scalar, 4>(powerSeriesCoeffs)
735 case unknownReactionRateType:
738 <<
"Internal error on line " << lineNo_-1
739 <<
": reaction rate type has not been set" 746 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
747 <<
" on line " << lineNo_-1
748 <<
" not implemented" 756 if (
mag(nAtoms[i]) > imbalanceTol_)
759 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
760 <<
" in " << elementNames_[i]
761 <<
" in reaction" <<
nl 762 << reactions_.last() <<
nl 763 <<
" on line " << lineNo_-1
770 reactionCoeffsTable.clear();
774 void Foam::chemkinReader::read
776 const fileName& CHEMKINFileName,
777 const fileName& thermoFileName,
778 const fileName& transportFileName
781 transportDict_.read(IFstream(transportFileName)());
785 std::ifstream thermoStream(thermoFileName.c_str());
790 <<
"file " << thermoFileName <<
" not found" 794 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
795 yy_switch_to_buffer(bufferPtr);
800 yy_delete_buffer(bufferPtr);
805 std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
810 <<
"file " << CHEMKINFileName <<
" not found" 814 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
815 yy_switch_to_buffer(bufferPtr);
817 initReactionKeywordTable();
822 yy_delete_buffer(bufferPtr);
828 Foam::chemkinReader::chemkinReader
839 speciesTable_(species),
840 reactions_(speciesTable_, speciesThermo_),
841 newFormat_(newFormat),
842 imbalanceTol_(ROOTSMALL)
844 read(CHEMKINFileName, thermoFileName, transportFileName);
848 Foam::chemkinReader::chemkinReader
856 speciesTable_(species),
857 reactions_(speciesTable_, speciesThermo_),
859 imbalanceTol_(thermoDict.
lookupOrDefault(
"imbalanceTolerance", ROOTSMALL))
863 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
870 if (thermoDict.
found(
"CHEMKINThermoFile"))
884 if (!chemkinFile.isAbsolute())
886 chemkinFile = relPath/chemkinFile;
891 thermoFile = relPath/thermoFile;
894 if (!transportFile.isAbsolute())
896 transportFile = relPath/transportFile;
900 read(chemkinFile, thermoFile, transportFile);
#define forAll(list, i)
Loop across all elements in list.
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.
const fileName & name() const
Return the dictionary name.
Macros for easy insertion into run-time selection tables.
bool read(const char *, int32_t &)
List< scalar > scalarList
A List of scalars.
string expand(const string &, const HashTable< string, word, string::hash > &mapping, const char sigil= '$')
Expand occurences of variables according to the mapping.
addChemistryReaderType(chemkinReader, gasHThermoPhysics)
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
atomicWeightTable atomicWeights
bool found(const Key &) const
Return true if hashedEntry is found in table.
bool isAbsolute() const
Return true if file name is absolute.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
A wordList with hashed indices for faster lookup by name.
dimensioned< scalar > mag(const dimensioned< Type > &)
const scalar RR
Universal gas constant (default in [J/(kmol K)])
fileName path() const
Return directory path name (part before last /)
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
sutherlandTransport< species::thermo< janafThermo< perfectGas< specie > >, sensibleEnthalpy > > gasHThermoPhysics
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.