59 "nonEquilibriumReversible",
67 "unimolecularFallOff",
68 "chemicallyActivatedBimolecular",
72 "unknownReactionRateType"
80 "unknownFallOffFunctionType"
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
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>
173 const ReactionRateType& rr
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>
249 const scalar Afactor0,
250 const scalar AfactorInf,
254 checkCoeffs(k0Coeffs,
"k0", 3);
255 checkCoeffs(kInfCoeffs,
"kInf", 3);
265 PressureDependencyType
270 Afactor0*k0Coeffs[0],
276 AfactorInf*kInfCoeffs[0],
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)
307 TroeCoeffs[3] = GREAT;
314 PressureDependencyType
319 Afactor0*k0Coeffs[0],
325 AfactorInf*kInfCoeffs[0],
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)
370 PressureDependencyType
375 Afactor0*k0Coeffs[0],
381 AfactorInf*kInfCoeffs[0],
401 <<
"Fall-off function type "
402 << fallOffFunctionNames[fofType]
403 <<
" on line " << lineNo_-1
404 <<
" not implemented"
424 checkCoeffs(ArrheniusCoeffs,
"Arrhenius", 3);
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;
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);
500 Afactor*ArrheniusCoeffs[0],
502 ArrheniusCoeffs[2]/
RR
506 AfactorRev*reverseArrheniusCoeffs[0],
507 reverseArrheniusCoeffs[1],
508 reverseArrheniusCoeffs[2]/
RR
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);
556 Afactor*concFactor*ArrheniusCoeffs[0],
558 ArrheniusCoeffs[2]/
RR,
563 AfactorRev*concFactor*reverseArrheniusCoeffs[0],
564 reverseArrheniusCoeffs[1],
565 reverseArrheniusCoeffs[2]/
RR,
579 Afactor*concFactor*ArrheniusCoeffs[0],
581 ArrheniusCoeffs[2]/
RR,
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);
658 Afactor*ArrheniusCoeffs[0],
660 ArrheniusCoeffs[2]/
RR,
661 LandauTellerCoeffs[0],
662 LandauTellerCoeffs[1]
666 AfactorRev*reverseArrheniusCoeffs[0],
667 reverseArrheniusCoeffs[1],
668 reverseArrheniusCoeffs[2]/
RR,
669 reverseLandauTellerCoeffs[0],
670 reverseLandauTellerCoeffs[1]
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,
717 reactionCoeffsTable[reactionRateTypeNames[rrType]];
719 checkCoeffs(powerSeriesCoeffs,
"power-series", 4);
727 Afactor*ArrheniusCoeffs[0],
729 ArrheniusCoeffs[2]/
RR,
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();
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);
839 speciesTable_(species),
840 reactions_(speciesTable_, speciesThermo_),
841 newFormat_(newFormat),
842 imbalanceTol_(ROOTSMALL)
844 read(CHEMKINFileName, thermoFileName, transportFileName);
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"))
886 chemkinFile = relPath/chemkinFile;
891 thermoFile = relPath/thermoFile;
896 transportFile = relPath/transportFile;
900 read(chemkinFile, thermoFile, transportFile);