DAC.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2016-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "DAC.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class ThermoType>
32 (
33  const IOdictionary& dict,
35 )
36 :
38  searchInitSet_(),
39  zprime_(0),
40  nbCLarge_(3),
41  sC_(this->nSpecie(),0),
42  sH_(this->nSpecie(),0),
43  sO_(this->nSpecie(),0),
44  sN_(this->nSpecie(),0),
45  CO2Id_(-1),
46  COId_(-1),
47  HO2Id_(-1),
48  H2OId_(-1),
49  NOId_(-1),
50  automaticSIS_(true),
51  phiTol_(this->tolerance()),
52  NOxThreshold_(1800),
53  CO2Name_
54  (
55  dict.subDict("reduction").lookupOrDefault<word>
56  (
57  "CO2Name","CO2"
58  )
59  ),
60  COName_
61  (
62  dict.subDict("reduction").lookupOrDefault<word>
63  (
64  "COName","CO"
65  )
66  ),
67  HO2Name_
68  (
69  dict.subDict("reduction").lookupOrDefault<word>
70  (
71  "HO2Name","HO2"
72  )
73  ),
74  H2OName_
75  (
76  dict.subDict("reduction").lookupOrDefault<word>
77  (
78  "H2OName","H2O"
79  )
80  ),
81  NOName_
82  (
83  dict.subDict("reduction").lookupOrDefault<word>
84  (
85  "NOName","NO"
86  )
87  ),
88  forceFuelInclusion_(false)
89 {
90  const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
91  forAllConstIter(wordHashSet, initSet, iter)
92  {
93  searchInitSet_.append(chemistry.mixture().species()[iter.key()]);
94  }
95 
96  if (this->coeffsDict_.found("automaticSIS"))
97  {
98  automaticSIS_.readIfPresent("automaticSIS", this->coeffsDict_);
99  }
100 
101  if (this->coeffsDict_.found("forceFuelInclusion"))
102  {
103  forceFuelInclusion_.readIfPresent
104  (
105  "forceFuelInclusion",this->coeffsDict_
106  );
107  }
108 
109  if (this->coeffsDict_.found("phiTol"))
110  {
111  phiTol_ = this->coeffsDict_.template lookup<scalar>("phiTol");
112  }
113 
114  if (this->coeffsDict_.found("NOxThreshold"))
115  {
116  NOxThreshold_ =
117  this->coeffsDict_.template lookup<scalar>("NOxThreshold");
118  }
119 
120  for (label i=0; i<this->nSpecie(); i++)
121  {
122  const List<specieElement>& curSpecieComposition =
123  chemistry.mixture().specieComposition(i);
124 
125  // For all elements in the current species
126  forAll(curSpecieComposition, j)
127  {
128  const specieElement& curElement =
129  curSpecieComposition[j];
130  if (curElement.name() == "C")
131  {
132  sC_[i] = curElement.nAtoms();
133  }
134  else if (curElement.name() == "H")
135  {
136  sH_[i] = curElement.nAtoms();
137  }
138  else if (curElement.name() == "O")
139  {
140  sO_[i] = curElement.nAtoms();
141  }
142  else if (curElement.name() == "N")
143  {
144  sN_[i] = curElement.nAtoms();
145  }
146  else
147  {
148  Info<< "element not considered"<<endl;
149  }
150  }
151  if (this->chemistry_.Y()[i].member() == CO2Name_)
152  {
153  CO2Id_ = i;
154  }
155  else if (this->chemistry_.Y()[i].member() == COName_)
156  {
157  COId_ = i;
158  }
159  else if (this->chemistry_.Y()[i].member() == HO2Name_)
160  {
161  HO2Id_ = i;
162  }
163  else if (this->chemistry_.Y()[i].member() == H2OName_)
164  {
165  H2OId_ = i;
166  }
167  else if (this->chemistry_.Y()[i].member() == NOName_)
168  {
169  NOId_ = i;
170  }
171  }
172 
173  if ((CO2Id_==-1 || COId_==-1 || HO2Id_==-1 || H2OId_==-1) && automaticSIS_)
174  {
176  << "The name of the species used in automatic SIS are not found in "
177  << " the mechanism. You should either set the name for CO2, CO, H2O"
178  << " and HO2 properly or set automaticSIS to off "
179  << exit(FatalError);
180  }
181 
182  // To compute zprime, the fuel species should be specified.
183  // According to the given mass fraction, an equivalent O/C ratio is computed
184  if (automaticSIS_)
185  {
186  List<Tuple2<word, scalar>> fuelSpeciesEntry
187  (
188  this->coeffsDict_.lookup("fuelSpecies")
189  );
190 
191  fuelSpecies_.setSize(fuelSpeciesEntry.size());
192  fuelSpeciesID_.setSize(fuelSpeciesEntry.size());
193  fuelSpeciesProp_.setSize(fuelSpeciesEntry.size());
194  scalar Mmtot(0.0);
195 
196  forAll(fuelSpeciesEntry, i)
197  {
198  fuelSpecies_[i] = fuelSpeciesEntry[i].first();
199  fuelSpeciesProp_[i] = fuelSpeciesEntry[i].second();
200  fuelSpeciesID_[i] =
201  this->chemistry_.mixture().species()[fuelSpecies_[i]];
202  scalar curMm =
203  this->chemistry_.specieThermos()[fuelSpeciesID_[i]].W();
204  Mmtot += fuelSpeciesProp_[i]/curMm;
205  }
206 
207  this->coeffsDict_.readIfPresent("nbCLarge", nbCLarge_);
208 
209  Mmtot = 1.0/Mmtot;
210  scalar nbC(0.0);
211  scalar nbO(0.0);
212  forAll(fuelSpecies_, i)
213  {
214  label curID = fuelSpeciesID_[i];
215  scalar curMm = this->chemistry_.specieThermos()[curID].W();
216 
217  nbC += fuelSpeciesProp_[i]*Mmtot/curMm*sC_[curID];
218  nbO += fuelSpeciesProp_[i]*Mmtot/curMm*sO_[curID];
219  }
220  zprime_ = nbO/nbC;
221  }
222 }
223 
224 
225 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
226 
227 template<class ThermoType>
229 {}
230 
231 
232 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
233 
234 template<class ThermoType>
236 (
237  const scalar p,
238  const scalar T,
239  const scalarField& c,
240  List<label>& ctos,
241  DynamicList<label>& stoc,
242  const label li
243 )
244 {
246 
247  scalarField c1(this->chemistry_.nEqns(), 0.0);
248  for(label i=0; i<this->nSpecie(); i++)
249  {
250  c1[i] = c[i];
251  }
252 
253  c1[this->nSpecie()] = T;
254  c1[this->nSpecie()+1] = p;
255 
256  // Compute the rAB matrix
257  RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
258  scalarField PA(this->nSpecie(),0.0);
259  scalarField CA(this->nSpecie(),0.0);
260 
261  // Number of initialised rAB for each lines
262  Field<label> NbrABInit(this->nSpecie(),0);
263  // Position of the initialised rAB, -1 when not initialised
264  RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
265  // Index of the other species involved in the rABNum
266  RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
267 
268  forAll(this->chemistry_.reactions(), i)
269  {
270  const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
271 
272  // for each reaction compute omegai
273  scalar omegaf, omegar;
274  const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
275 
276  // Then for each pair of species composing this reaction,
277  // compute the rAB matrix (separate the numerator and
278  // denominator)
279 
280  // While computing the rAB for all the species involved in the reaction
281  // we should consider that one can write a reaction A+B=2C as A+B=C+C
282  // In this case, the following algorithm only take once the effect
283  // of the species. It stores the species encountered in the reaction but
284  // use another list to see if this species has already been used
285 
286  DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
287  DynamicList<label> wAID(R.lhs().size()+R.rhs().size());
288 
289  forAll(R.lhs(), s) // Compute rAB for all species in the left hand side
290  {
291  label ss = R.lhs()[s].index;
292  scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
293  List<bool> deltaBi(this->nSpecie(), false);
294  FIFOStack<label> usedIndex;
295  forAll(R.lhs(), j)
296  {
297  label sj = R.lhs()[j].index;
298  usedIndex.push(sj);
299  deltaBi[sj] = true;
300  }
301  forAll(R.rhs(), j)
302  {
303  label sj = R.rhs()[j].index;
304  usedIndex.push(sj);
305  deltaBi[sj] = true;
306  }
307 
308  // Disable for self reference (by definition rAA=0)
309  deltaBi[ss] = false;
310 
311  while(!usedIndex.empty())
312  {
313  label curIndex = usedIndex.pop();
314  if (deltaBi[curIndex])
315  {
316  // Disable to avoid counting it more than once
317  deltaBi[curIndex] = false;
318  // Test if this rAB is not initialised
319  if (rABPos(ss, curIndex)==-1)
320  {
321  // It starts at rABPos(ss, sj)=0
322  rABPos(ss, curIndex) = NbrABInit[ss];
323  NbrABInit[ss]++;
324  // to avoid overflow
325  rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
326  // store the other specie involved
327  rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
328  }
329  else
330  {
331  rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
332  }
333  }
334  }
335 
336  bool found(false);
337  forAll(wAID, id)
338  {
339  if (ss==wAID[id])
340  {
341  wA[id] += sl*omegai;
342  found = true;
343  }
344  }
345  if (!found)
346  {
347  wA.append(sl*omegai);
348  wAID.append(ss);
349  }
350  }
351 
352  forAll(R.rhs(), s) // Compute rAB for all species in the right hand side
353  {
354  label ss = R.rhs()[s].index;
355  scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
356  List<bool> deltaBi(this->nSpecie(), false);
357  FIFOStack<label> usedIndex;
358  forAll(R.lhs(), j)
359  {
360  label sj = R.lhs()[j].index;
361  usedIndex.push(sj);
362  deltaBi[sj] = true;
363  }
364  forAll(R.rhs(), j)
365  {
366  label sj = R.rhs()[j].index;
367  usedIndex.push(sj);
368  deltaBi[sj] = true;
369  }
370 
371  // Disable for self reference (by definition rAA=0)
372  deltaBi[ss] = false;
373 
374  while(!usedIndex.empty())
375  {
376  label curIndex = usedIndex.pop();
377  if (deltaBi[curIndex])
378  {
379  // Disable to avoid counting it more than once
380  deltaBi[curIndex] = false;
381 
382  // Test if this rAB is not initialised
383  if (rABPos(ss, curIndex) == -1)
384  {
385  // it starts at rABPos(ss, sj)=0
386  rABPos(ss, curIndex) = NbrABInit[ss];
387  NbrABInit[ss]++;
388  rABNum(ss, rABPos(ss, curIndex)) = sl*omegai;
389  rABOtherSpec(ss, rABPos(ss, curIndex)) = curIndex;
390  }
391  else
392  {
393  rABNum(ss, rABPos(ss, curIndex)) += sl*omegai;
394  }
395  }
396  }
397 
398  bool found(false);
399  forAll(wAID, id)
400  {
401  if (ss==wAID[id])
402  {
403  wA[id] += sl*omegai;
404  found = true;
405  }
406  }
407  if (!found)
408  {
409  wA.append(sl*omegai);
410  wAID.append(ss);
411  }
412  }
413  wAID.shrink();
414 
415  // Now that every species of the reactions has been visited, we can
416  // compute the production and consumption rate. This way, it avoids
417  // getting wrong results when species are present in both lhs and rhs
418  forAll(wAID, id)
419  {
420  if (wA[id] > 0.0)
421  {
422  if (PA[wAID[id]] == 0.0)
423  {
424  PA[wAID[id]] = wA[id];
425  }
426  else
427  {
428  PA[wAID[id]] += wA[id];
429  }
430  }
431  else
432  {
433  if (CA[wAID[id]] == 0.0)
434  {
435  CA[wAID[id]] = -wA[id];
436  }
437  else
438  {
439  CA[wAID[id]] += -wA[id];
440  }
441  }
442  }
443  }
444  // rii = 0.0 by definition
445 
446  scalar phiLarge(0.0);
447  scalar phiProgress(0.0);
448  if (automaticSIS_)
449  {
450  // Compute the progress equivalence ratio
451  // and the equivalence ratio for fuel decomposition
452  label nElements = 4; // 4 main elements (C, H, O, N)
453 
454  // Total number of C, H and O (in this order)
455  scalarList Na(nElements,0.0);
456  scalarList Nal(nElements,0.0); // for large hydrocarbons
457 
458  for (label i=0; i<this->nSpecie(); i++)
459  {
460  // Complete combustion products are not considered
461  if
462  (
463  this->chemistry_.Y()[i].member() == "CO2"
464  || this->chemistry_.Y()[i].member() == "H2O"
465  )
466  {
467  continue;
468  }
469  Na[0] += sC_[i]*c[i];
470  Na[1] += sH_[i]*c[i];
471  Na[2] += sO_[i]*c[i];
472  if (sC_[i]>nbCLarge_ || this->chemistry_.Y()[i].member() == "O2")
473  {
474  Nal[0] += sC_[i]*c[i];
475  Nal[1] += sH_[i]*c[i];
476  Nal[2] += sO_[i]*c[i];
477  }
478  }
479 
480  // 2C(-CO2) + H(-H2O)/2 - z'C(-CO2)
481  // Progress equivalence ratio = ----------------------------------
482  // O(-CO2-H2O) - z' C(-CO2)
483  // where minus means that this species is not considered for the number
484  // of atoms and z' is the ratio of the number of O and C in the fuel(s)
485  phiProgress = (2*Na[0]+Na[1]/2-zprime_*Na[0])/(Na[2]-zprime_*Na[0]);
486 
487  // 2Cl + Hl/2
488  // Equivalence ratio for fuel decomposition = ----------
489  // Ol(+O2)
490  phiLarge = (2*Nal[0]+Nal[1]/2)/Nal[2];
491  }
492 
493  // Using the rAB matrix (numerator and denominator separated)
494  // compute the R value according to the search initiating set
495  scalarField Rvalue(this->nSpecie(),0.0);
496 
497  // Set all species to inactive and activate them according
498  // to rAB and initial set
499  for (label i=0; i<this->nSpecie(); i++)
500  {
501  this->activeSpecies_[i] = false;
502  }
503 
504  // Initialise the FIFOStack for search set
506 
507  const labelList& SIS(searchInitSet_);
508 
509  // If automaticSIS is on, the search initiating set is selected according to
510  // phiProgress and phiLarge
511  if (automaticSIS_)
512  {
513  if (phiLarge >= phiTol_ && phiProgress >= phiTol_)
514  {
515  // When phiLarge and phiProgress >= phiTol then
516  // CO, HO2 and fuel are in the SIS
517  Q.push(COId_);
518  this->activeSpecies_[COId_] = true;
519  Rvalue[COId_] = 1.0;
520  Q.push(HO2Id_);
521  this->activeSpecies_[HO2Id_] = true;
522  Rvalue[HO2Id_] = 1.0;
523  forAll(fuelSpeciesID_,i)
524  {
525  Q.push(fuelSpeciesID_[i]);
526  this->activeSpecies_[fuelSpeciesID_[i]] = true;
527  Rvalue[fuelSpeciesID_[i]] = 1.0;
528  }
529 
530  }
531  else if (phiLarge < phiTol_ && phiProgress >= phiTol_)
532  {
533  // When phiLarge < phiTol and phiProgress >= phiTol then
534  // CO, HO2 are in the SIS
535  Q.push(COId_);
536  this->activeSpecies_[COId_] = true;
537  Rvalue[COId_] = 1.0;
538  Q.push(HO2Id_);
539  this->activeSpecies_[HO2Id_] = true;
540  Rvalue[HO2Id_] = 1.0;
541 
542  if (forceFuelInclusion_)
543  {
544  forAll(fuelSpeciesID_,i)
545  {
546  Q.push(fuelSpeciesID_[i]);
547  this->activeSpecies_[fuelSpeciesID_[i]] = true;
548  Rvalue[fuelSpeciesID_[i]] = 1.0;
549  }
550  }
551  }
552  else
553  {
554  // When phiLarge and phiProgress< phiTol then
555  // CO2, H2O are in the SIS
556  Q.push(CO2Id_);
557  this->activeSpecies_[CO2Id_] = true;
558  Rvalue[CO2Id_] = 1.0;
559 
560  Q.push(H2OId_);
561  this->activeSpecies_[H2OId_] = true;
562  Rvalue[H2OId_] = 1.0;
563  if (forceFuelInclusion_)
564  {
565  forAll(fuelSpeciesID_,i)
566  {
567  Q.push(fuelSpeciesID_[i]);
568  this->activeSpecies_[fuelSpeciesID_[i]] = true;
569  Rvalue[fuelSpeciesID_[i]] = 1.0;
570  }
571  }
572  }
573 
574  if (T>NOxThreshold_ && NOId_!=-1)
575  {
576  Q.push(NOId_);
577  this->activeSpecies_[NOId_] = true;
578  Rvalue[NOId_] = 1.0;
579  }
580  }
581  else // No automaticSIS => all species of the SIS are added
582  {
583  for (label i=0; i<SIS.size(); i++)
584  {
585  label q = SIS[i];
586  this->activeSpecies_[q] = true;
587  Q.push(q);
588  Rvalue[q] = 1.0;
589  }
590  }
591 
592  // Execute the main loop for R-value
593  while (!Q.empty())
594  {
595  label u = Q.pop();
596  scalar Den = max(PA[u],CA[u]);
597  if (Den != 0)
598  {
599  for (label v=0; v<NbrABInit[u]; v++)
600  {
601  label otherSpec = rABOtherSpec(u, v);
602  scalar rAB = mag(rABNum(u, v))/Den;
603  if (rAB > 1)
604  {
605  rAB = 1;
606  }
607  // The direct link is weaker than the user-defined tolerance
608  if (rAB >= this->tolerance())
609  {
610  scalar Rtemp = Rvalue[u]*rAB;
611  // a link analysed previously is stronger
612  // the (composed) link is stronger than the user-defined
613  // tolerance
614  if ((Rvalue[otherSpec]<Rtemp) && (Rtemp>=this->tolerance()))
615  {
616  Q.push(otherSpec);
617  Rvalue[otherSpec] = Rtemp;
618  if (!this->activeSpecies_[otherSpec])
619  {
620  this->activeSpecies_[otherSpec] = true;
621  }
622  }
623  }
624  }
625  }
626  }
627 
629 }
630 
631 
632 // ************************************************************************* //
A FIFO stack based on a singly-linked list.
Definition: FIFOStack.H:51
A HashTable with keys but without contents.
Definition: HashSet.H:59
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
scalar omega(const scalar p, const scalar T, const scalarField &c, const label li, scalar &omegaf, scalar &omegar) const
Net reaction rate.
Definition: Reaction.C:319
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const label nSpecie
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: reactionI.H:36
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:55
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
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))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
A class for handling words, derived from string.
Definition: word.H:59
const word & name() const
Return the name of the element.
label nAtoms() const
Return the number of atoms of this element in the specie.
basicChemistryModel & chemistry
const multiComponentMixture< ThermoType > & mixture() const
Return reference to the mixture.
virtual ~DAC()
Destructor.
Definition: DAC.C:228
void setSize(const label)
Reset size of List.
Definition: List.C:281
const volScalarField & T
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: reactionI.H:42
#define R(A, B, C, D, E, F, K, M)
DAC(const IOdictionary &dict, chemistryModel< ThermoType > &chemistry)
Construct from components.
Definition: DAC.C:32
void push(const T &a)
Push an element onto the stack.
Definition: FIFOStack.H:84
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
T pop()
Pop the bottom element off the stack.
Definition: FIFOStack.H:90
An abstract class for methods of chemical mechanism reduction.
volScalarField & p
bool found
virtual void reduceMechanism(const scalar p, const scalar T, const scalarField &c, List< label > &ctos, DynamicList< label > &stoc, const label li)
Reduce the mechanism.
Definition: DAC.C:236