32 template<
class CompType,
class ThermoType>
40 searchInitSet_(this->coeffsDict_.subDict(
"initialSet").size()),
42 nbCLarge_(this->coeffsDict_.template lookupOrDefault<label>(
"nbCLarge", 3)),
43 sC_(this->nSpecie_, 0),
44 sH_(this->nSpecie_, 0),
45 sO_(this->nSpecie_, 0),
46 sN_(this->nSpecie_, 0),
54 this->coeffsDict_.template lookupOrDefault<Switch>
62 this->coeffsDict_.template lookupOrDefault<scalar>
64 "phiTol", this->tolerance()
69 this->coeffsDict_.template lookupOrDefault<scalar>
75 CO2Name_(this->coeffsDict_.template lookupOrDefault<word>(
"CO2",
"CO2")),
76 COName_(this->coeffsDict_.template lookupOrDefault<word>(
"CO",
"CO")),
77 HO2Name_(this->coeffsDict_.template lookupOrDefault<word>(
"HO2",
"HO2")),
78 H2OName_(this->coeffsDict_.template lookupOrDefault<word>(
"H2O",
"H2O")),
79 NOName_(this->coeffsDict_.template lookupOrDefault<word>(
"NO",
"NO")),
82 this->coeffsDict_.template lookupOrDefault<Switch>
96 searchInitSet_[j++] = i;
100 if (j < searchInitSet_.size())
103 << searchInitSet_.size()-j
104 <<
" species in the intial set is not in the mechanism "
113 for (
label i=0; i<this->nSpecie_; i++)
116 specieComposition[i];
119 forAll(curSpecieComposition, j)
123 if (curElement.
name() ==
"C")
125 sC_[i] = curElement.
nAtoms();
127 else if (curElement.
name() ==
"H")
129 sH_[i] = curElement.
nAtoms();
131 else if (curElement.
name() ==
"O")
133 sO_[i] = curElement.
nAtoms();
135 else if (curElement.
name() ==
"N")
137 sN_[i] = curElement.
nAtoms();
141 Info<<
" element " << curElement.
name() <<
" not considered"
145 if (this->chemistry_.Y()[i].member() == CO2Name_)
149 else if (this->chemistry_.Y()[i].member() == COName_)
153 else if (this->chemistry_.Y()[i].member() == HO2Name_)
157 else if (this->chemistry_.Y()[i].member() == H2OName_)
161 else if (this->chemistry_.Y()[i].member() == NOName_)
167 if ((CO2Id_==-1 || COId_==-1 || HO2Id_==-1 || H2OId_==-1) && automaticSIS_)
170 <<
"The name of the species used in automatic SIS are not found in "
171 <<
" the mechanism. You should either set the name for CO2, CO, H2O"
172 <<
" and HO2 properly or set automaticSIS to off "
181 if (this->coeffsDict_.found(
"fuelSpecies"))
183 fuelDict = this->coeffsDict_.
subDict(
"fuelSpecies");
184 fuelSpecies_ = fuelDict.
toc();
185 if (fuelSpecies_.size() == 0)
188 <<
"With automatic SIS, the fuel species should be "
189 <<
"specified in the fuelSpecies subDict"
196 <<
"With automatic SIS, the fuel species should be "
197 <<
"specified in the fuelSpecies subDict"
202 fuelSpeciesID_.setSize(fuelSpecies_.size());
203 fuelSpeciesProp_.setSize(fuelSpecies_.size());
208 fuelSpeciesProp_[i] = fuelDict.
get<scalar>(fuelSpecies_[i]);
209 for (
label j=0; j<this->nSpecie_; j++)
211 if (this->chemistry_.Y()[j].member() == fuelSpecies_[i])
213 fuelSpeciesID_[i] = j;
218 this->chemistry_.specieThermo()[fuelSpeciesID_[i]].W();
219 Mmtot += fuelSpeciesProp_[i]/curMm;
227 label curID = fuelSpeciesID_[i];
228 scalar curMm = this->chemistry_.specieThermo()[curID].W();
230 nbC += fuelSpeciesProp_[i]*Mmtot/curMm*sC_[curID];
231 nbO += fuelSpeciesProp_[i]*Mmtot/curMm*sO_[curID];
240 template<
class CompType,
class ThermoType>
248 template<
class CompType,
class ThermoType>
256 scalarField& completeC(this->chemistry_.completeC());
258 for (
label i=0; i<this->nSpecie_; i++)
264 c1[this->nSpecie_] =
T;
265 c1[this->nSpecie_+1] =
p;
279 scalar pf, cf, pr, cr;
284 scalar omegai = this->chemistry_.omega
286 R,
c1,
T,
p, pf, cf, lRef, pr, cr, rRef
305 scalar sl = -
R.lhs()[
s].stoichCoeff;
310 label sj =
R.lhs()[j].index;
316 label sj =
R.rhs()[j].index;
324 while (!usedIndex.empty())
327 if (deltaBi[curIndex])
330 deltaBi[curIndex] =
false;
332 if (rABPos(ss, curIndex) == -1)
335 rABPos(ss, curIndex) = NbrABInit[ss];
338 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
340 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
344 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
360 wA.append(sl*omegai);
369 scalar sl =
R.rhs()[
s].stoichCoeff;
374 label sj =
R.lhs()[j].index;
380 label sj =
R.rhs()[j].index;
388 while (!usedIndex.empty())
391 if (deltaBi[curIndex])
394 deltaBi[curIndex] =
false;
397 if (rABPos(ss, curIndex) == -1)
400 rABPos(ss, curIndex) = NbrABInit[ss];
402 rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
403 rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
407 rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
423 wA.append(sl*omegai);
436 if (PA[wAID[
id]] == 0.0)
438 PA[wAID[id]] = wA[id];
442 PA[wAID[id]] += wA[id];
447 if (CA[wAID[
id]] == 0.0)
449 CA[wAID[id]] = -wA[id];
453 CA[wAID[id]] += -wA[id];
460 scalar phiLarge(0.0);
461 scalar phiProgress(0.0);
472 for (
label i=0; i<this->nSpecie_; i++)
477 this->chemistry_.Y()[i].member() ==
"CO2"
478 || this->chemistry_.Y()[i].member() ==
"H2O"
484 Na[0] += sC_[i]*
c[i];
485 Na[1] += sH_[i]*
c[i];
486 Na[2] += sO_[i]*
c[i];
487 if (sC_[i]>nbCLarge_ || this->chemistry_.Y()[i].member() ==
"O2")
489 Nal[0] += sC_[i]*
c[i];
490 Nal[1] += sH_[i]*
c[i];
491 Nal[2] += sO_[i]*
c[i];
500 phiProgress = (2*Na[0]+Na[1]/2-zprime_*Na[0])/(Na[2]-zprime_*Na[0]);
505 phiLarge = (2*Nal[0] + Nal[1]/2)/Nal[2];
511 label speciesNumber = 0;
515 for (
label i=0; i<this->nSpecie_; ++i)
517 this->activeSpecies_[i] =
false;
529 if (phiLarge >= phiTol_ && phiProgress >= phiTol_)
535 this->activeSpecies_[COId_] =
true;
539 this->activeSpecies_[HO2Id_] =
true;
540 Rvalue[HO2Id_] = 1.0;
541 for (
const label fuelId : fuelSpeciesID_)
545 this->activeSpecies_[fuelId] =
true;
546 Rvalue[fuelId] = 1.0;
550 else if (phiLarge < phiTol_ && phiProgress >= phiTol_)
556 this->activeSpecies_[COId_] =
true;
560 this->activeSpecies_[HO2Id_] =
true;
561 Rvalue[HO2Id_] = 1.0;
563 if (forceFuelInclusion_)
565 for (
const label fuelId : fuelSpeciesID_)
569 this->activeSpecies_[fuelId] =
true;
570 Rvalue[fuelId] = 1.0;
580 this->activeSpecies_[CO2Id_] =
true;
581 Rvalue[CO2Id_] = 1.0;
585 this->activeSpecies_[H2OId_] =
true;
586 Rvalue[H2OId_] = 1.0;
587 if (forceFuelInclusion_)
589 for (
const label fuelId : fuelSpeciesID_)
593 this->activeSpecies_[fuelId] =
true;
594 Rvalue[fuelId] = 1.0;
599 if (
T > NOxThreshold_ && NOId_ != -1)
603 this->activeSpecies_[NOId_] =
true;
609 for (
label i=0; i<SIS.size(); i++)
612 this->activeSpecies_[q] =
true;
623 scalar Den =
max(PA[u], CA[u]);
626 for (
label v=0; v<NbrABInit[u]; ++v)
628 label otherSpec = rABOtherSpec(u, v);
629 scalar rAB =
mag(rABNum(u, v))/Den;
636 if (rAB >= this->tolerance())
638 scalar Rtemp = Rvalue[u]*rAB;
642 if ((Rvalue[otherSpec]<Rtemp) && (Rtemp>=this->tolerance()))
645 Rvalue[otherSpec] = Rtemp;
646 if (!this->activeSpecies_[otherSpec])
648 this->activeSpecies_[otherSpec] =
true;
658 forAll(this->chemistry_.reactions(), i)
661 this->chemistry_.reactionsDisabled()[i] =
false;
666 if (!this->activeSpecies_[ss])
669 this->chemistry_.reactionsDisabled()[i] =
true;
673 if (!this->chemistry_.reactionsDisabled()[i])
678 if (!this->activeSpecies_[ss])
681 this->chemistry_.reactionsDisabled()[i] =
true;
688 this->NsSimp_ = speciesNumber;
689 scalarField& simplifiedC(this->chemistry_.simplifiedC());
690 simplifiedC.setSize(this->NsSimp_ + 2);
693 Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
696 for (
label i=0; i<this->nSpecie_; i++)
698 if (this->activeSpecies_[i])
701 simplifiedC[j] =
c[i];
703 if (!this->chemistry_.active(i))
705 this->chemistry_.setActive(i);
714 simplifiedC[this->NsSimp_] =
T;
715 simplifiedC[this->NsSimp_+1] =
p;
716 this->chemistry_.setNsDAC(this->NsSimp_);
720 this->chemistry_.setNSpecie(this->NsSimp_);