58 const char* Foam::chemkinReader::reactionTypeNames[4] =
62 "nonEquilibriumReversible",
66 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
70 "unimolecularFallOff",
71 "chemicallyActivatedBimolecular",
75 "unknownReactionRateType" 78 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
83 "unknownFallOffFunctionType" 86 void Foam::chemkinReader::initReactionKeywordTable()
88 reactionKeywordTable_.
insert(
"M", thirdBodyReactionType);
89 reactionKeywordTable_.
insert(
"LOW", unimolecularFallOffReactionType);
90 reactionKeywordTable_.
insert 93 chemicallyActivatedBimolecularReactionType
95 reactionKeywordTable_.
insert(
"TROE", TroeReactionType);
96 reactionKeywordTable_.
insert(
"SRI", SRIReactionType);
97 reactionKeywordTable_.
insert(
"LT", LandauTellerReactionType);
98 reactionKeywordTable_.
insert(
"RLT", reverseLandauTellerReactionType);
99 reactionKeywordTable_.
insert(
"JAN", JanevReactionType);
100 reactionKeywordTable_.
insert(
"FIT1", powerSeriesReactionRateType);
101 reactionKeywordTable_.
insert(
"HV", radiationActivatedReactionType);
102 reactionKeywordTable_.
insert(
"TDEP", speciesTempReactionType);
103 reactionKeywordTable_.
insert(
"EXCI", energyLossReactionType);
104 reactionKeywordTable_.
insert(
"MOME", plasmaMomentumTransfer);
105 reactionKeywordTable_.
insert(
"XSMI", collisionCrossSection);
106 reactionKeywordTable_.
insert(
"REV", nonEquilibriumReversibleReactionType);
107 reactionKeywordTable_.
insert(
"DUPLICATE", duplicateReactionType);
108 reactionKeywordTable_.
insert(
"DUP", duplicateReactionType);
109 reactionKeywordTable_.
insert(
"FORD", speciesOrderForward);
110 reactionKeywordTable_.
insert(
"RORD", speciesOrderReverse);
111 reactionKeywordTable_.
insert(
"UNITS", UnitsOfReaction);
112 reactionKeywordTable_.
insert(
"END", end);
116 Foam::scalar Foam::chemkinReader::molecularWeight
118 const List<specieElement>& specieComposition
123 forAll(specieComposition, i)
125 label nAtoms = specieComposition[i].nAtoms();
126 const word& elementName = specieComposition[i].name();
128 if (isotopeAtomicWts_.found(elementName))
130 molWt += nAtoms*isotopeAtomicWts_[elementName];
139 <<
"Unknown element " << elementName
140 <<
" on line " << lineNo_-1 <<
nl 141 <<
" specieComposition: " << specieComposition
150 void Foam::chemkinReader::checkCoeffs
153 const char* reactionRateName,
157 if (reactionCoeffs.size() != nCoeffs)
160 <<
"Wrong number of coefficients for the " << reactionRateName
161 <<
" rate expression on line " 162 << lineNo_-1 <<
", should be " 163 << nCoeffs <<
" but " << reactionCoeffs.size() <<
" supplied." <<
nl 164 <<
"Coefficients are " 165 << reactionCoeffs <<
nl 170 template<
class ReactionRateType>
171 void Foam::chemkinReader::addReactionType
173 const reactionType rType,
174 DynamicList<gasHReaction::specieCoeffs>& lhs,
175 DynamicList<gasHReaction::specieCoeffs>& rhs,
176 const ReactionRateType& rr
185 new IrreversibleReaction
188 Reaction<gasHThermoPhysics>
205 new ReversibleReaction
208 Reaction<gasHThermoPhysics>
226 <<
"Reaction type " << reactionTypeNames[rType]
227 <<
" on line " << lineNo_-1
228 <<
" not handled by this function" 234 <<
"Unknown reaction type " << rType
235 <<
" on line " << lineNo_-1
241 template<
template<
class,
class>
class PressureDependencyType>
242 void Foam::chemkinReader::addPressureDependentReaction
244 const reactionType rType,
245 const fallOffFunctionType fofType,
246 DynamicList<gasHReaction::specieCoeffs>& lhs,
247 DynamicList<gasHReaction::specieCoeffs>& rhs,
251 const HashTable<scalarList>& reactionCoeffsTable,
252 const scalar Afactor0,
253 const scalar AfactorInf,
257 checkCoeffs(k0Coeffs,
"k0", 3);
258 checkCoeffs(kInfCoeffs,
"kInf", 3);
268 PressureDependencyType
269 <ArrheniusReactionRate, LindemannFallOffFunction>
271 ArrheniusReactionRate
273 Afactor0*k0Coeffs[0],
277 ArrheniusReactionRate
279 AfactorInf*kInfCoeffs[0],
283 LindemannFallOffFunction(),
284 thirdBodyEfficiencies(speciesTable_, efficiencies)
293 reactionCoeffsTable[fallOffFunctionNames[fofType]]
296 if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
299 <<
"Wrong number of coefficients for Troe rate expression" 300 " on line " << lineNo_-1 <<
", should be 3 or 4 but " 301 << TroeCoeffs.size() <<
" supplied." <<
nl 302 <<
"Coefficients are " 307 if (TroeCoeffs.size() == 3)
309 TroeCoeffs.setSize(4);
310 TroeCoeffs[3] = GREAT;
317 PressureDependencyType
318 <ArrheniusReactionRate, TroeFallOffFunction>
320 ArrheniusReactionRate
322 Afactor0*k0Coeffs[0],
326 ArrheniusReactionRate
328 AfactorInf*kInfCoeffs[0],
339 thirdBodyEfficiencies(speciesTable_, efficiencies)
348 reactionCoeffsTable[fallOffFunctionNames[fofType]]
351 if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
354 <<
"Wrong number of coefficients for SRI rate expression" 355 " on line " << lineNo_-1 <<
", should be 3 or 5 but " 356 << SRICoeffs.size() <<
" supplied." <<
nl 357 <<
"Coefficients are " 362 if (SRICoeffs.size() == 3)
364 SRICoeffs.setSize(5);
373 PressureDependencyType
374 <ArrheniusReactionRate, SRIFallOffFunction>
376 ArrheniusReactionRate
378 Afactor0*k0Coeffs[0],
382 ArrheniusReactionRate
384 AfactorInf*kInfCoeffs[0],
396 thirdBodyEfficiencies(speciesTable_, efficiencies)
404 <<
"Fall-off function type " 405 << fallOffFunctionNames[fofType]
406 <<
" on line " << lineNo_-1
407 <<
" not implemented" 414 void Foam::chemkinReader::addReaction
416 DynamicList<gasHReaction::specieCoeffs>& lhs,
417 DynamicList<gasHReaction::specieCoeffs>& rhs,
419 const reactionType rType,
420 const reactionRateType rrType,
421 const fallOffFunctionType fofType,
423 HashTable<scalarList>& reactionCoeffsTable,
427 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
433 const List<specieElement>& specieComposition =
434 speciesComposition_[speciesTable_[lhs[i].index]];
436 forAll(specieComposition, j)
438 label elementi = elementIndices_[specieComposition[j].name()];
440 lhs[i].stoichCoeff*specieComposition[j].nAtoms();
446 const List<specieElement>& specieComposition =
447 speciesComposition_[speciesTable_[rhs[i].index]];
449 forAll(specieComposition, j)
451 label elementi = elementIndices_[specieComposition[j].name()];
453 rhs[i].stoichCoeff*specieComposition[j].nAtoms();
460 const scalar concFactor = 0.001;
464 sumExp += lhs[i].exponent;
466 scalar Afactor =
pow(concFactor, sumExp - 1.0);
468 scalar AfactorRev = Afactor;
470 if (rType == nonEquilibriumReversible)
475 sumExp += rhs[i].exponent;
477 AfactorRev =
pow(concFactor, sumExp - 1.0);
484 if (rType == nonEquilibriumReversible)
487 reactionCoeffsTable[reactionTypeNames[rType]];
489 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
493 new NonEquilibriumReversibleReaction
496 Reaction<gasHThermoPhysics>
503 ArrheniusReactionRate
505 Afactor*ArrheniusCoeffs[0],
507 ArrheniusCoeffs[2]/
RR 509 ArrheniusReactionRate
511 AfactorRev*reverseArrheniusCoeffs[0],
512 reverseArrheniusCoeffs[1],
513 reverseArrheniusCoeffs[2]/
RR 524 ArrheniusReactionRate
526 Afactor*ArrheniusCoeffs[0],
528 ArrheniusCoeffs[2]/
RR 534 case thirdBodyArrhenius:
536 if (rType == nonEquilibriumReversible)
539 reactionCoeffsTable[reactionTypeNames[rType]];
541 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
545 new NonEquilibriumReversibleReaction
549 thirdBodyArrheniusReactionRate
552 Reaction<gasHThermoPhysics>
559 thirdBodyArrheniusReactionRate
561 Afactor*concFactor*ArrheniusCoeffs[0],
563 ArrheniusCoeffs[2]/
RR,
564 thirdBodyEfficiencies(speciesTable_, efficiencies)
566 thirdBodyArrheniusReactionRate
568 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
569 reverseArrheniusCoeffs[1],
570 reverseArrheniusCoeffs[2]/
RR,
571 thirdBodyEfficiencies(speciesTable_, efficiencies)
582 thirdBodyArrheniusReactionRate
584 Afactor*concFactor*ArrheniusCoeffs[0],
586 ArrheniusCoeffs[2]/
RR,
587 thirdBodyEfficiencies(speciesTable_, efficiencies)
593 case unimolecularFallOff:
595 addPressureDependentReaction<FallOffReactionRate>
602 reactionCoeffsTable[reactionRateTypeNames[rrType]],
611 case chemicallyActivatedBimolecular:
613 addPressureDependentReaction<ChemicallyActivatedReactionRate>
621 reactionCoeffsTable[reactionRateTypeNames[rrType]],
632 reactionCoeffsTable[reactionRateTypeNames[rrType]];
633 checkCoeffs(LandauTellerCoeffs,
"Landau-Teller", 2);
635 if (rType == nonEquilibriumReversible)
638 reactionCoeffsTable[reactionTypeNames[rType]];
639 checkCoeffs(reverseArrheniusCoeffs,
"reverse Arrhenius", 3);
644 word(reactionTypeNames[rType])
645 + reactionRateTypeNames[rrType]
647 checkCoeffs(LandauTellerCoeffs,
"reverse Landau-Teller", 2);
651 new NonEquilibriumReversibleReaction
654 Reaction<gasHThermoPhysics>
661 LandauTellerReactionRate
663 Afactor*ArrheniusCoeffs[0],
665 ArrheniusCoeffs[2]/
RR,
666 LandauTellerCoeffs[0],
667 LandauTellerCoeffs[1]
669 LandauTellerReactionRate
671 AfactorRev*reverseArrheniusCoeffs[0],
672 reverseArrheniusCoeffs[1],
673 reverseArrheniusCoeffs[2]/
RR,
674 reverseLandauTellerCoeffs[0],
675 reverseLandauTellerCoeffs[1]
686 LandauTellerReactionRate
688 Afactor*ArrheniusCoeffs[0],
690 ArrheniusCoeffs[2]/
RR,
691 LandauTellerCoeffs[0],
692 LandauTellerCoeffs[1]
701 reactionCoeffsTable[reactionRateTypeNames[rrType]];
703 checkCoeffs(JanevCoeffs,
"Janev", 9);
711 Afactor*ArrheniusCoeffs[0],
713 ArrheniusCoeffs[2]/
RR,
714 FixedList<scalar, 9>(JanevCoeffs)
722 reactionCoeffsTable[reactionRateTypeNames[rrType]];
724 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
730 powerSeriesReactionRate
732 Afactor*ArrheniusCoeffs[0],
734 ArrheniusCoeffs[2]/
RR,
735 FixedList<scalar, 4>(powerSeriesCoeffs)
740 case unknownReactionRateType:
743 <<
"Internal error on line " << lineNo_-1
744 <<
": reaction rate type has not been set" 751 <<
"Reaction rate type " << reactionRateTypeNames[rrType]
752 <<
" on line " << lineNo_-1
753 <<
" not implemented" 761 if (
mag(nAtoms[i]) > imbalanceTol_)
764 <<
"Elemental imbalance of " <<
mag(nAtoms[i])
765 <<
" in " << elementNames_[i]
766 <<
" in reaction" <<
nl 767 << reactions_.last() <<
nl 768 <<
" on line " << lineNo_-1
775 reactionCoeffsTable.clear();
779 void Foam::chemkinReader::read
781 const fileName& CHEMKINFileName,
782 const fileName& thermoFileName,
783 const fileName& transportFileName
786 transportDict_.read(IFstream(transportFileName)());
788 if (!thermoFileName.empty())
790 std::ifstream thermoStream(thermoFileName);
795 <<
"file " << thermoFileName <<
" not found" 799 yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
800 yy_switch_to_buffer(bufferPtr);
805 yy_delete_buffer(bufferPtr);
810 std::ifstream CHEMKINStream(CHEMKINFileName);
815 <<
"file " << CHEMKINFileName <<
" not found" 819 yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
820 yy_switch_to_buffer(bufferPtr);
822 initReactionKeywordTable();
827 yy_delete_buffer(bufferPtr);
833 Foam::chemkinReader::chemkinReader
844 speciesTable_(species),
845 reactions_(speciesTable_, speciesThermo_),
846 newFormat_(newFormat),
847 imbalanceTol_(ROOTSMALL)
849 read(CHEMKINFileName, thermoFileName, transportFileName);
853 Foam::chemkinReader::chemkinReader
861 speciesTable_(species),
862 reactions_(speciesTable_, speciesThermo_),
863 newFormat_(
thermoDict.getOrDefault(
"newFormat", false)),
864 imbalanceTol_(
thermoDict.getOrDefault(
"imbalanceTolerance", ROOTSMALL))
868 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
871 fileName chemkinFile(
thermoDict.get<fileName>(
"CHEMKINFile"));
872 chemkinFile.expand();
875 thermoDict.readIfPresent(
"CHEMKINThermoFile", thermoFile);
878 fileName transportFile(
thermoDict.get<fileName>(
"CHEMKINTransportFile"));
879 transportFile.expand();
885 if (!chemkinFile.isAbsolute())
887 chemkinFile = relPath/chemkinFile;
890 if (!thermoFile.empty() && !thermoFile.isAbsolute())
892 thermoFile = relPath/thermoFile;
895 if (!transportFile.isAbsolute())
897 transportFile = relPath/transportFile;
901 read(chemkinFile, thermoFile, transportFile);
List< scalar > scalarList
List of scalar.
A class for handling file names.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
bool found(const Key &key) const
Same as contains()
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Macros for easy insertion into run-time selection tables.
bool read(const char *buf, int32_t &val)
Same as readInt32.
#define forAll(list, i)
Loop across all elements in list.
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
addChemistryReaderType(chemkinReader, gasHThermoPhysics)
const dictionary & thermoDict
atomicWeightTable atomicWeights
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
A wordList with hashed named lookup, which can be faster in some situations than using the normal lis...
messageStream Info
Information stream (stdout output on master, null elsewhere)
const scalar RR
Universal gas constant: default in [J/(kmol K)].
sutherlandTransport< species::thermo< janafThermo< perfectGas< specie > >, sensibleEnthalpy > > gasHThermoPhysics
static constexpr const zero Zero
Global zero (0)