30 template<
class CompType,
class ThermoType>
38 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size()),
41 sC_(this->nSpecie_,0),
42 sH_(this->nSpecie_,0),
43 sO_(this->nSpecie_,0),
44 sN_(this->nSpecie_,0),
51 phiTol_(this->tolerance()),
55 dict.subDict(
"reduction").lookupOrDefault<
word>
62 dict.subDict(
"reduction").lookupOrDefault<
word>
69 dict.subDict(
"reduction").lookupOrDefault<
word>
76 dict.subDict(
"reduction").lookupOrDefault<
word>
83 dict.subDict(
"reduction").lookupOrDefault<
word>
88 forceFuelInclusion_(
false)
94 if (initSet.
found(chemistry.
Y()[i].member()))
96 searchInitSet_[j++] = i;
99 if (j<searchInitSet_.size())
102 << searchInitSet_.size()-j
103 <<
" species in the initial set is not in the mechanism " 108 if (this->coeffsDict_.found(
"automaticSIS"))
110 automaticSIS_.readIfPresent(
"automaticSIS", this->coeffsDict_);
113 if (this->coeffsDict_.found(
"forceFuelInclusion"))
115 forceFuelInclusion_.readIfPresent
117 "forceFuelInclusion",this->coeffsDict_
121 if (this->coeffsDict_.found(
"phiTol"))
123 phiTol_ =
readScalar(this->coeffsDict_.lookup(
"phiTol"));
126 if (this->coeffsDict_.found(
"NOxThreshold"))
128 NOxThreshold_ =
readScalar(this->coeffsDict_.lookup(
"NOxThreshold"));
133 for (
label i=0; i<this->nSpecie_; i++)
136 specieComposition[i];
138 forAll(curSpecieComposition, j)
141 curSpecieComposition[j];
142 if (curElement.
name() ==
"C")
144 sC_[i] = curElement.
nAtoms();
146 else if (curElement.
name() ==
"H")
148 sH_[i] = curElement.
nAtoms();
150 else if (curElement.
name() ==
"O")
152 sO_[i] = curElement.
nAtoms();
154 else if (curElement.
name() ==
"N")
156 sN_[i] = curElement.
nAtoms();
160 Info<<
"element not considered"<<
endl;
163 if (this->chemistry_.Y()[i].member() == CO2Name_)
167 else if (this->chemistry_.Y()[i].member() == COName_)
171 else if (this->chemistry_.Y()[i].member() == HO2Name_)
175 else if (this->chemistry_.Y()[i].member() == H2OName_)
179 else if (this->chemistry_.Y()[i].member() == NOName_)
185 if ((CO2Id_==-1 || COId_==-1 || HO2Id_==-1 || H2OId_==-1) && automaticSIS_)
188 <<
"The name of the species used in automatic SIS are not found in " 189 <<
" the mechanism. You should either set the name for CO2, CO, H2O" 190 <<
" and HO2 properly or set automaticSIS to off " 199 if (this->coeffsDict_.found(
"fuelSpecies"))
201 fuelDict = this->coeffsDict_.
subDict(
"fuelSpecies");
202 fuelSpecies_ = fuelDict.
toc();
203 if (fuelSpecies_.size() == 0)
206 <<
"With automatic SIS, the fuel species should be " 207 <<
"specified in the fuelSpecies subDict" 214 <<
"With automatic SIS, the fuel species should be " 215 <<
"specified in the fuelSpecies subDict" 219 if (this->coeffsDict_.found(
"nbCLarge"))
224 fuelSpeciesID_.setSize(fuelSpecies_.size());
225 fuelSpeciesProp_.setSize(fuelSpecies_.size());
231 for (
label j=0; j<this->nSpecie_; j++)
233 if (this->chemistry_.Y()[j].member() == fuelSpecies_[i])
235 fuelSpeciesID_[i] = j;
240 this->chemistry_.specieThermo()[fuelSpeciesID_[i]].W();
241 Mmtot += fuelSpeciesProp_[i]/curMm;
249 label curID = fuelSpeciesID_[i];
250 scalar curMm = this->chemistry_.specieThermo()[curID].W();
252 nbC += fuelSpeciesProp_[i]*Mmtot/curMm*sC_[curID];
253 nbO += fuelSpeciesProp_[i]*Mmtot/curMm*sO_[curID];
262 template<
class CompType,
class ThermoType>
270 template<
class CompType,
class ThermoType>
278 scalarField& completeC(this->chemistry_.completeC());
280 for(
label i=0; i<this->nSpecie_; i++)
286 c1[this->nSpecie_] =
T;
287 c1[this->nSpecie_+1] =
p;
301 scalar pf, cf, pr, cr;
303 forAll(this->chemistry_.reactions(), i)
307 scalar omegai = this->chemistry_.omega
309 R, c1, T, p, pf, cf, lRef, pr, cr, rRef
328 scalar sl = -R.
lhs()[
s].stoichCoeff;
347 while(!usedIndex.empty())
350 if (deltaBi[curIndex])
353 deltaBi[curIndex] =
false;
355 if (rABPos(ss, curIndex)==-1)
358 rABPos(ss, curIndex) = NbrABInit[ss];
361 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
363 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
367 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
383 wA.append(sl*omegai);
391 scalar sl = R.
rhs()[
s].stoichCoeff;
410 while(!usedIndex.empty())
413 if (deltaBi[curIndex])
416 deltaBi[curIndex] =
false;
419 if (rABPos(ss, curIndex) == -1)
422 rABPos(ss, curIndex) = NbrABInit[ss];
424 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
425 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
429 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
445 wA.append(sl*omegai);
458 if (PA[wAID[
id]] == 0.0)
460 PA[wAID[id]] = wA[id];
464 PA[wAID[id]] += wA[id];
469 if (CA[wAID[
id]] == 0.0)
471 CA[wAID[id]] = -wA[id];
475 CA[wAID[id]] += -wA[id];
482 scalar phiLarge(0.0);
483 scalar phiProgress(0.0);
494 for (
label i=0; i<this->nSpecie_; i++)
499 this->chemistry_.Y()[i].member() ==
"CO2" 500 || this->chemistry_.Y()[i].member() ==
"H2O" 505 Na[0] += sC_[i]*c[i];
506 Na[1] += sH_[i]*c[i];
507 Na[2] += sO_[i]*c[i];
508 if (sC_[i]>nbCLarge_ || this->chemistry_.Y()[i].member() ==
"O2")
510 Nal[0] += sC_[i]*c[i];
511 Nal[1] += sH_[i]*c[i];
512 Nal[2] += sO_[i]*c[i];
521 phiProgress = (2*Na[0]+Na[1]/2-zprime_*Na[0])/(Na[2]-zprime_*Na[0]);
526 phiLarge = (2*Nal[0]+Nal[1]/2)/Nal[2];
532 label speciesNumber = 0;
536 for (
label i=0; i<this->nSpecie_; i++)
538 this->activeSpecies_[i] =
false;
550 if (phiLarge >= phiTol_ && phiProgress >= phiTol_)
556 this->activeSpecies_[COId_] =
true;
560 this->activeSpecies_[HO2Id_] =
true;
561 Rvalue[HO2Id_] = 1.0;
564 Q.
push(fuelSpeciesID_[i]);
566 this->activeSpecies_[fuelSpeciesID_[i]] =
true;
567 Rvalue[fuelSpeciesID_[i]] = 1.0;
571 else if (phiLarge < phiTol_ && phiProgress >= phiTol_)
577 this->activeSpecies_[COId_] =
true;
581 this->activeSpecies_[HO2Id_] =
true;
582 Rvalue[HO2Id_] = 1.0;
584 if (forceFuelInclusion_)
588 Q.
push(fuelSpeciesID_[i]);
590 this->activeSpecies_[fuelSpeciesID_[i]] =
true;
591 Rvalue[fuelSpeciesID_[i]] = 1.0;
601 this->activeSpecies_[CO2Id_] =
true;
602 Rvalue[CO2Id_] = 1.0;
606 this->activeSpecies_[H2OId_] =
true;
607 Rvalue[H2OId_] = 1.0;
608 if (forceFuelInclusion_)
612 Q.
push(fuelSpeciesID_[i]);
614 this->activeSpecies_[fuelSpeciesID_[i]] =
true;
615 Rvalue[fuelSpeciesID_[i]] = 1.0;
620 if (T>NOxThreshold_ && NOId_!=-1)
624 this->activeSpecies_[NOId_] =
true;
633 this->activeSpecies_[q] =
true;
644 scalar Den =
max(PA[u],CA[u]);
647 for (
label v=0; v<NbrABInit[u]; v++)
649 label otherSpec = rABOtherSpec(u, v);
650 scalar rAB =
mag(rABNum(u, v))/Den;
656 if (rAB >= this->tolerance())
658 scalar Rtemp = Rvalue[u]*rAB;
662 if ((Rvalue[otherSpec]<Rtemp) && (Rtemp>=this->tolerance()))
665 Rvalue[otherSpec] = Rtemp;
666 if (!this->activeSpecies_[otherSpec])
668 this->activeSpecies_[otherSpec] =
true;
678 forAll(this->chemistry_.reactions(), i)
681 this->chemistry_.reactionsDisabled()[i] =
false;
686 if (!this->activeSpecies_[ss])
689 this->chemistry_.reactionsDisabled()[i] =
true;
693 if (!this->chemistry_.reactionsDisabled()[i])
698 if (!this->activeSpecies_[ss])
701 this->chemistry_.reactionsDisabled()[i] =
true;
708 this->NsSimp_ = speciesNumber;
709 scalarField& simplifiedC(this->chemistry_.simplifiedC());
710 simplifiedC.
setSize(this->NsSimp_+2);
713 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
716 for (
label i=0; i<this->nSpecie_; i++)
718 if (this->activeSpecies_[i])
721 simplifiedC[j] = c[i];
723 if (!this->chemistry_.active(i))
725 this->chemistry_.setActive(i);
734 simplifiedC[this->NsSimp_] =
T;
735 simplifiedC[this->NsSimp_+1] =
p;
736 this->chemistry_.setNsDAC(this->NsSimp_);
740 this->chemistry_.setNSpecie(this->NsSimp_);
A FIFO stack based on a singly-linked list.
virtual label nSpecie() const
The number of species.
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.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
PtrList< volScalarField > & Y()
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...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void size(const label)
Override size to be inconsistent with allocated storage.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
BasicChemistryModel< rhoReactionThermo > & chemistry
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
wordList toc() const
Return the table of contents.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
void setSize(const label)
Alter the addressed list size.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
A class for handling words, derived from string.
DAC(const IOdictionary &dict, TDACChemistryModel< CompType, ThermoType > &chemistry)
Construct from components.
const word & name() const
Return the name of the element.
label nAtoms() const
Return the number of atoms of this element in the specie.
label readLabel(Istream &is)
List< List< specieElement > > & specieComp()
virtual ~DAC()
Destructor.
virtual void reduceMechanism(const scalarField &c, const scalar T, const scalar p)
Reduce the mechanism.
void setSize(const label)
Reset size of List.
#define R(A, B, C, D, E, F, K, M)
void push(const T &a)
Push an element onto the stack.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
dimensioned< scalar > mag(const dimensioned< Type > &)
T pop()
Pop the bottom element off the stack.
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
An abstract class for methods of chemical mechanism reduction.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.