pyrolysisChemistryModel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2013-2017 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 
27 #include "solidReaction.H"
28 #include "basicThermo.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class CompType, class SolidThermo, class GasThermo>
35 (
36  const fvMesh& mesh,
37  const word& phaseName
38 )
39 :
41  pyrolisisGases_(this->reactions_[0].gasSpecies()),
42  gasThermo_(pyrolisisGases_.size()),
43  nGases_(pyrolisisGases_.size()),
44  nSpecie_(this->Ys_.size() + nGases_),
45  RRg_(nGases_),
46  Ys0_(this->nSolids_),
47  cellCounter_(0)
48 {
49  // create the fields for the chemistry sources
50  forAll(this->RRs_, fieldi)
51  {
52  IOobject header
53  (
54  this->Ys_[fieldi].name() + "0",
55  mesh.time().timeName(),
56  mesh,
57  IOobject::NO_READ
58  );
59 
60  // check if field exists and can be read
61  if (header.typeHeaderOk<volScalarField>(true))
62  {
63  Ys0_.set
64  (
65  fieldi,
66  new volScalarField
67  (
68  IOobject
69  (
70  this->Ys_[fieldi].name() + "0",
71  mesh.time().timeName(),
72  mesh,
73  IOobject::MUST_READ,
74  IOobject::AUTO_WRITE
75  ),
76  mesh
77  )
78  );
79  }
80  else
81  {
82  volScalarField Y0Default
83  (
84  IOobject
85  (
86  "Y0Default",
87  mesh.time().timeName(),
88  mesh,
89  IOobject::MUST_READ,
90  IOobject::NO_WRITE
91  ),
92  mesh
93  );
94 
95  Ys0_.set
96  (
97  fieldi,
98  new volScalarField
99  (
100  IOobject
101  (
102  this->Ys_[fieldi].name() + "0",
103  mesh.time().timeName(),
104  mesh,
105  IOobject::NO_READ,
106  IOobject::AUTO_WRITE
107  ),
108  Y0Default
109  )
110  );
111 
112  // Calculate inital values of Ysi0 = rho*delta*Yi
113  Ys0_[fieldi].primitiveFieldRef() =
114  this->solidThermo().rho()
115  *max(this->Ys_[fieldi], scalar(0.001))*mesh.V();
116  }
117  }
118 
119  forAll(RRg_, fieldi)
120  {
121  RRg_.set
122  (
123  fieldi,
125  (
126  IOobject
127  (
128  "RRg." + pyrolisisGases_[fieldi],
129  mesh.time().timeName(),
130  mesh,
131  IOobject::NO_READ,
132  IOobject::NO_WRITE
133  ),
134  mesh,
136  )
137  );
138  }
139 
140  forAll(gasThermo_, gasI)
141  {
142  dictionary thermoDict =
143  mesh.lookupObject<dictionary>
144  (
146  ).subDict(pyrolisisGases_[gasI]);
147 
148  gasThermo_.set
149  (
150  gasI,
151  new GasThermo(thermoDict)
152  );
153  }
154 
155  Info<< "pyrolysisChemistryModel: " << nl;
156  Info<< indent << "Number of solids = " << this->nSolids_ << nl;
157  Info<< indent << "Number of gases = " << nGases_ << nl;
158  forAll(this->reactions_, i)
159  {
160  Info<< dynamic_cast<const solidReaction<SolidThermo>& >
161  (
162  this->reactions_[i]
163  ) << nl;
164  }
165 }
166 
167 
168 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
169 
170 template<class CompType, class SolidThermo, class GasThermo>
173 {}
174 
175 
176 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
177 
178 template<class CompType, class SolidThermo, class GasThermo>
181 (
182  const scalarField& c,
183  const scalar T,
184  const scalar p,
185  const bool updateC0
186 ) const
187 {
188  scalar pf, cf, pr, cr;
189  label lRef, rRef;
190 
191  const label celli = cellCounter_;
192 
193  scalarField om(nEqns(), 0.0);
194 
195  forAll(this->reactions_, i)
196  {
197  const Reaction<SolidThermo>& R = this->reactions_[i];
198 
199  scalar omegai = omega
200  (
201  R, c, T, p, pf, cf, lRef, pr, cr, rRef
202  );
203  scalar rhoL = 0.0;
204  forAll(R.lhs(), s)
205  {
206  label si = R.lhs()[s].index;
207  om[si] -= omegai;
208  rhoL = this->solidThermo_[si].rho(p, T);
209  }
210  scalar sr = 0.0;
211  forAll(R.rhs(), s)
212  {
213  label si = R.rhs()[s].index;
214  scalar rhoR = this->solidThermo_[si].rho(p, T);
215  sr = rhoR/rhoL;
216  om[si] += sr*omegai;
217 
218  if (updateC0)
219  {
220  Ys0_[si][celli] += sr*omegai;
221  }
222  }
223  forAll(R.grhs(), g)
224  {
225  label gi = R.grhs()[g].index;
226  om[gi + this->nSolids_] += (1.0 - sr)*omegai;
227  }
228  }
229 
230  return om;
231 }
232 
233 
234 template<class CompType, class SolidThermo, class GasThermo>
235 Foam::scalar
237 (
238  const Reaction<SolidThermo>& R,
239  const scalarField& c,
240  const scalar T,
241  const scalar p,
242  scalar& pf,
243  scalar& cf,
244  label& lRef,
245  scalar& pr,
246  scalar& cr,
247  label& rRef
248 ) const
249 {
250  scalarField c1(nSpecie_, 0.0);
251 
252  label celli = cellCounter_;
253 
254  for (label i=0; i<nSpecie_; i++)
255  {
256  c1[i] = max(0.0, c[i]);
257  }
258 
259  scalar kf = R.kf(p, T, c1);
260 
261  const label Nl = R.lhs().size();
262 
263  for (label s=0; s<Nl; s++)
264  {
265  label si = R.lhs()[s].index;
266  const scalar exp = R.lhs()[si].exponent;
267 
268  kf *=
269  pow(c1[si]/Ys0_[si][celli], exp)
270  *(Ys0_[si][celli]);
271  }
272 
273  return kf;
274 }
275 
276 
277 template<class CompType, class SolidThermo, class GasThermo>
279 omegaI
280 (
281  const label index,
282  const scalarField& c,
283  const scalar T,
284  const scalar p,
285  scalar& pf,
286  scalar& cf,
287  label& lRef,
288  scalar& pr,
289  scalar& cr,
290  label& rRef
291 ) const
292 {
293 
294  const Reaction<SolidThermo>& R = this->reactions_[index];
295  scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
296  return(w);
297 }
298 
299 
300 template<class CompType, class SolidThermo, class GasThermo>
303 (
304  const scalar time,
305  const scalarField &c,
306  scalarField& dcdt
307 ) const
308 {
309  const scalar T = c[nSpecie_];
310  const scalar p = c[nSpecie_ + 1];
311 
312  dcdt = 0.0;
313 
314  dcdt = omega(c, T, p);
315 
316  //Total mass concentration
317  scalar cTot = 0.0;
318  for (label i=0; i<this->nSolids_; i++)
319  {
320  cTot += c[i];
321  }
322 
323  scalar newCp = 0.0;
324  scalar newhi = 0.0;
325  for (label i=0; i<this->nSolids_; i++)
326  {
327  scalar dYidt = dcdt[i]/cTot;
328  scalar Yi = c[i]/cTot;
329  newCp += Yi*this->solidThermo_[i].Cp(p, T);
330  newhi -= dYidt*this->solidThermo_[i].Hc();
331  }
332 
333  scalar dTdt = newhi/newCp;
334  scalar dtMag = min(500.0, mag(dTdt));
335  dcdt[nSpecie_] = dTdt*dtMag/(mag(dTdt) + 1.0e-10);
336 
337  // dp/dt = ...
338  dcdt[nSpecie_ + 1] = 0.0;
339 }
340 
341 
342 template<class CompType, class SolidThermo, class GasThermo>
344 jacobian
345 (
346  const scalar t,
347  const scalarField& c,
348  scalarField& dcdt,
349  scalarSquareMatrix& dfdc
350 ) const
351 {
352  const scalar T = c[nSpecie_];
353  const scalar p = c[nSpecie_ + 1];
354 
355  scalarField c2(nSpecie_, 0.0);
356 
357  for (label i=0; i<this->nSolids_; i++)
358  {
359  c2[i] = max(c[i], 0.0);
360  }
361 
362  for (label i=0; i<nEqns(); i++)
363  {
364  for (label j=0; j<nEqns(); j++)
365  {
366  dfdc(i, j) = 0.0;
367  }
368  }
369 
370  // length of the first argument must be nSolids
371  dcdt = omega(c2, T, p);
372 
373  for (label ri=0; ri<this->reactions_.size(); ri++)
374  {
375  const Reaction<SolidThermo>& R = this->reactions_[ri];
376 
377  scalar kf0 = R.kf(p, T, c2);
378 
379  forAll(R.lhs(), j)
380  {
381  label sj = R.lhs()[j].index;
382  scalar kf = kf0;
383  forAll(R.lhs(), i)
384  {
385  label si = R.lhs()[i].index;
386  scalar exp = R.lhs()[i].exponent;
387  if (i == j)
388  {
389  if (exp < 1.0)
390  {
391  if (c2[si]>SMALL)
392  {
393  kf *= exp*pow(c2[si] + VSMALL, exp - 1.0);
394  }
395  else
396  {
397  kf = 0.0;
398  }
399  }
400  else
401  {
402  kf *= exp*pow(c2[si], exp - 1.0);
403  }
404  }
405  else
406  {
407  Info<< "Solid reactions have only elements on slhs"
408  << endl;
409  kf = 0.0;
410  }
411  }
412 
413  forAll(R.lhs(), i)
414  {
415  label si = R.lhs()[i].index;
416  dfdc[si][sj] -= kf;
417  }
418  forAll(R.rhs(), i)
419  {
420  label si = R.rhs()[i].index;
421  dfdc[si][sj] += kf;
422  }
423  }
424  }
425 
426  // calculate the dcdT elements numerically
427  scalar delta = 1.0e-8;
428  scalarField dcdT0 = omega(c2, T - delta, p);
429  scalarField dcdT1 = omega(c2, T + delta, p);
430 
431  for (label i=0; i<nEqns(); i++)
432  {
433  dfdc[i][nSpecie_] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
434  }
435 
436 }
437 
438 
439 template<class CompType, class SolidThermo, class GasThermo>
442 {
443  // nEqns = number of solids + gases + temperature + pressure
444  return (nSpecie_ + 2);
445 }
446 
447 
448 template<class CompType, class SolidThermo, class GasThermo>
451 {
452  if (!this->chemistry_)
453  {
454  return;
455  }
456 
457  const volScalarField rho
458  (
459  IOobject
460  (
461  "rho",
462  this->time().timeName(),
463  this->mesh(),
464  IOobject::NO_READ,
465  IOobject::NO_WRITE,
466  false
467  ),
468  this->solidThermo().rho()
469  );
470 
471  forAll(this->RRs_, i)
472  {
473  this->RRs_[i].field() = 0.0;
474  }
475 
476  forAll(RRg_, i)
477  {
478  RRg_[i].field() = 0.0;
479  }
480 
481  forAll(rho, celli)
482  {
483  cellCounter_ = celli;
484 
485  const scalar delta = this->mesh().V()[celli];
486 
487  if (this->reactingCells_[celli])
488  {
489  scalar rhoi = rho[celli];
490  scalar Ti = this->solidThermo().T()[celli];
491  scalar pi = this->solidThermo().p()[celli];
492 
493  scalarField c(nSpecie_, 0.0);
494  for (label i=0; i<this->nSolids_; i++)
495  {
496  c[i] = rhoi*this->Ys_[i][celli]*delta;
497  }
498 
499  const scalarField dcdt = omega(c, Ti, pi, true);
500 
501  forAll(this->RRs_, i)
502  {
503  this->RRs_[i][celli] = dcdt[i]/delta;
504  }
505 
506  forAll(RRg_, i)
507  {
508  RRg_[i][celli] = dcdt[this->nSolids_ + i]/delta;
509  }
510  }
511  }
512 }
513 
514 
515 template<class CompType, class SolidThermo, class GasThermo>
516 Foam::scalar
518 (
519  const scalar deltaT
520 )
521 {
522  scalar deltaTMin = GREAT;
523 
524  if (!this->chemistry_)
525  {
526  return deltaTMin;
527  }
528 
529  const volScalarField rho
530  (
531  IOobject
532  (
533  "rho",
534  this->time().timeName(),
535  this->mesh(),
536  IOobject::NO_READ,
537  IOobject::NO_WRITE,
538  false
539  ),
540  this->solidThermo().rho()
541  );
542 
543  forAll(this->RRs_, i)
544  {
545  this->RRs_[i].field() = 0.0;
546  }
547  forAll(RRg_, i)
548  {
549  RRg_[i].field() = 0.0;
550  }
551 
552  const scalarField& T = this->solidThermo().T();
553  const scalarField& p = this->solidThermo().p();
554 
555  scalarField c(nSpecie_, 0.0);
556  scalarField c0(nSpecie_, 0.0);
557  scalarField dc(nSpecie_, 0.0);
558  scalarField delta(this->mesh().V());
559 
560  forAll(rho, celli)
561  {
562  if (this->reactingCells_[celli])
563  {
564  cellCounter_ = celli;
565 
566  scalar rhoi = rho[celli];
567  scalar pi = p[celli];
568  scalar Ti = T[celli];
569 
570  for (label i=0; i<this->nSolids_; i++)
571  {
572  c[i] = rhoi*this->Ys_[i][celli]*delta[celli];
573  }
574 
575  c0 = c;
576 
577  // Initialise time progress
578  scalar timeLeft = deltaT;
579 
580  // calculate the chemical source terms
581  while (timeLeft > SMALL)
582  {
583  scalar dt = timeLeft;
584  this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
585  timeLeft -= dt;
586  }
587 
588  deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
589  dc = c - c0;
590 
591  forAll(this->RRs_, i)
592  {
593  this->RRs_[i][celli] = dc[i]/(deltaT*delta[celli]);
594  }
595 
596  forAll(RRg_, i)
597  {
598  RRg_[i][celli] = dc[this->nSolids_ + i]/(deltaT*delta[celli]);
599  }
600 
601  // Update Ys0_
602  dc = omega(c0, Ti, pi, true);
603  }
604  }
605 
606  // Don't allow the time-step to change more than a factor of 2
607  deltaTMin = min(deltaTMin, 2*deltaT);
608 
609  return deltaTMin;
610 }
611 
612 
613 template<class CompType, class SolidThermo,class GasThermo>
616 (
617  const volScalarField& p,
618  const volScalarField& T,
619  const label index
620 ) const
621 {
623  (
624  new volScalarField
625  (
626  IOobject
627  (
628  "Hs_" + pyrolisisGases_[index],
629  this->mesh_.time().timeName(),
630  this->mesh_,
631  IOobject::NO_READ,
632  IOobject::NO_WRITE,
633  false
634  ),
635  this->mesh_,
636  dimensionedScalar("zero", dimEnergy/dimMass, 0.0)
637  )
638  );
639 
640  volScalarField::Internal& gasHs = tHs.ref();
641 
642  const GasThermo& mixture = gasThermo_[index];
643 
644  forAll(gasHs, celli)
645  {
646  gasHs[celli] = mixture.Hs(p[celli], T[celli]);
647  }
648 
649  return tHs;
650 }
651 
652 
653 template<class CompType, class SolidThermo, class GasThermo>
655 (
656  scalarField &c,
657  scalar& T,
658  scalar& p,
659  scalar& deltaT,
660  scalar& subDeltaT
661 ) const
662 {
664 }
665 
666 
667 // ************************************************************************* //
scalar delta
virtual scalar kf(const scalar p, const scalar T, const scalarField &c) const
Forward rate constant.
Definition: Reaction.C:392
virtual ~pyrolysisChemistryModel()
Destructor.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
const List< specieCoeffs > & lhs() const
Definition: ReactionI.H:51
virtual const List< specieCoeffs > & grhs() const
Definition: Reaction.C:453
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual scalar omegaI(label iReaction, const scalarField &c, const scalar T, const scalar p, scalar &pf, scalar &cf, label &lRef, scalar &pr, scalar &cr, label &rRef) const
Return the reaction rate for iReaction.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
Definition: solidThermo.C:120
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
pyrolysisChemistryModel(const fvMesh &mesh, const word &phaseName)
Construct from mesh and phase name.
virtual const volScalarField & T() const
Temperature [K].
Definition: basicThermo.C:489
Extends base solid chemistry model by adding a thermo package, and ODE functions. Introduces chemistr...
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:53
virtual volScalarField & p()
Pressure [Pa].
Definition: basicThermo.C:477
virtual scalarField omega(const scalarField &c, const scalar T, const scalar p, const bool updateC0=false) const
dc/dt = omega, rate of change in concentration, for each species
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
virtual label nEqns() const
Number of ODE&#39;s to solve.
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))
dynamicFvMesh & mesh
dimensionedScalar exp(const dimensionedScalar &ds)
virtual tmp< volScalarField > gasHs(const volScalarField &p, const volScalarField &T, const label i) const
Return sensible enthalpy for gas i [J/Kg].
A class for handling words, derived from string.
Definition: word.H:59
virtual scalar solve(const scalar deltaT)
Solve the reaction system for the given time step.
rhoEqn solve()
word timeName
Definition: getTimeIndex.H:3
const word dictName("particleTrackDict")
static const char nl
Definition: Ostream.H:262
const dimensionedVector & g
Fundamental solid thermodynamic properties.
Definition: solidThermo.H:54
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const dimensionSet dimEnergy
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const dimensionedScalar c
Speed of light in a vacuum.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
A class for managing temporary objects.
Definition: PtrList.H:53
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
virtual void calculate()
Calculates the reaction rates.
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
const List< specieCoeffs > & rhs() const
Definition: ReactionI.H:59