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 <<
"Unhandled reaction type " << reactionTypeNames[rType]
227 <<
" on line " << lineNo_-1 <<
nl 233 <<
"Unknown reaction type (" << int(rType)
234 <<
"), on line " << lineNo_-1 <<
nl 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 (" << int(fofType)
405 <<
") not implemented, line " << lineNo_-1
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 Reaction<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 Reaction<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
652 Reaction<gasHThermoPhysics>
659 LandauTellerReactionRate
661 Afactor*ArrheniusCoeffs[0],
663 ArrheniusCoeffs[2]/
RR,
664 LandauTellerCoeffs[0],
665 LandauTellerCoeffs[1]
667 LandauTellerReactionRate
669 AfactorRev*reverseArrheniusCoeffs[0],
670 reverseArrheniusCoeffs[1],
671 reverseArrheniusCoeffs[2]/
RR,
672 reverseLandauTellerCoeffs[0],
673 reverseLandauTellerCoeffs[1]
684 LandauTellerReactionRate
686 Afactor*ArrheniusCoeffs[0],
688 ArrheniusCoeffs[2]/
RR,
689 LandauTellerCoeffs[0],
690 LandauTellerCoeffs[1]
699 reactionCoeffsTable[reactionRateTypeNames[rrType]];
701 checkCoeffs(JanevCoeffs,
"Janev", 9);
709 Afactor*ArrheniusCoeffs[0],
711 ArrheniusCoeffs[2]/
RR,
712 FixedList<scalar, 9>(JanevCoeffs)
720 reactionCoeffsTable[reactionRateTypeNames[rrType]];
722 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
728 powerSeriesReactionRate
730 Afactor*ArrheniusCoeffs[0],
732 ArrheniusCoeffs[2]/
RR,
733 FixedList<scalar, 4>(powerSeriesCoeffs)
738 case unknownReactionRateType:
741 <<
"Internal error on line " << lineNo_-1
742 <<
": reaction rate type has not been set" 749 <<
"Reaction rate type (" << int(rrType)
750 <<
") not implemented, on line " << lineNo_-1
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)());
785 if (!thermoFileName.empty())
787 std::ifstream thermoStream(thermoFileName);
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);
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_),
860 newFormat_(
thermoDict.getOrDefault(
"newFormat", false)),
861 imbalanceTol_(
thermoDict.getOrDefault(
"imbalanceTolerance", ROOTSMALL))
865 Info<<
"Reading CHEMKIN thermo data in new file format" <<
endl;
868 fileName chemkinFile(
thermoDict.get<fileName>(
"CHEMKINFile"));
869 chemkinFile.expand();
872 thermoDict.readIfPresent(
"CHEMKINThermoFile", thermoFile);
875 fileName transportFile(
thermoDict.get<fileName>(
"CHEMKINTransportFile"));
876 transportFile.expand();
882 if (!chemkinFile.isAbsolute())
884 chemkinFile = relPath/chemkinFile;
887 if (!thermoFile.empty() && !thermoFile.isAbsolute())
889 thermoFile = relPath/thermoFile;
892 if (!transportFile.isAbsolute())
894 transportFile = relPath/transportFile;
898 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)