chemistryModel.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) 2011-2015 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 "chemistryModel.H"
27 #include "reactingMixture.H"
28 #include "UniformField.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class CompType, class ThermoType>
34 (
35  const fvMesh& mesh,
36  const word& phaseName
37 )
38 :
39  CompType(mesh, phaseName),
40  ODESystem(),
41  Y_(this->thermo().composition().Y()),
42  reactions_
43  (
44  dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
45  ),
46  specieThermo_
47  (
48  dynamic_cast<const reactingMixture<ThermoType>&>
49  (this->thermo()).speciesData()
50  ),
51 
52  nSpecie_(Y_.size()),
53  nReaction_(reactions_.size()),
54  Treact_(CompType::template lookupOrDefault<scalar>("Treact", 0.0)),
55  RR_(nSpecie_)
56 {
57  // create the fields for the chemistry sources
58  forAll(RR_, fieldI)
59  {
60  RR_.set
61  (
62  fieldI,
64  (
65  IOobject
66  (
67  "RR." + Y_[fieldI].name(),
68  mesh.time().timeName(),
69  mesh,
70  IOobject::NO_READ,
71  IOobject::NO_WRITE
72  ),
73  mesh,
75  )
76  );
77  }
78 
79  Info<< "chemistryModel: Number of species = " << nSpecie_
80  << " and reactions = " << nReaction_ << endl;
81 }
82 
83 
84 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
85 
86 template<class CompType, class ThermoType>
88 {}
89 
90 
91 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
92 
93 template<class CompType, class ThermoType>
96 (
97  const scalarField& c,
98  const scalar T,
99  const scalar p
100 ) const
101 {
102  scalar pf, cf, pr, cr;
103  label lRef, rRef;
104 
105  tmp<scalarField> tom(new scalarField(nEqns(), 0.0));
106  scalarField& om = tom();
107 
108  forAll(reactions_, i)
109  {
110  const Reaction<ThermoType>& R = reactions_[i];
111 
112  scalar omegai = omega
113  (
114  R, c, T, p, pf, cf, lRef, pr, cr, rRef
115  );
116 
117  forAll(R.lhs(), s)
118  {
119  const label si = R.lhs()[s].index;
120  const scalar sl = R.lhs()[s].stoichCoeff;
121  om[si] -= sl*omegai;
122  }
123 
124  forAll(R.rhs(), s)
125  {
126  const label si = R.rhs()[s].index;
127  const scalar sr = R.rhs()[s].stoichCoeff;
128  om[si] += sr*omegai;
129  }
130  }
131 
132  return tom;
133 }
134 
135 
136 template<class CompType, class ThermoType>
138 (
139  const label index,
140  const scalarField& c,
141  const scalar T,
142  const scalar p,
143  scalar& pf,
144  scalar& cf,
145  label& lRef,
146  scalar& pr,
147  scalar& cr,
148  label& rRef
149 ) const
150 {
151 
152  const Reaction<ThermoType>& R = reactions_[index];
153  scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
154  return(w);
155 }
156 
157 
158 template<class CompType, class ThermoType>
160 (
161  const Reaction<ThermoType>& R,
162  const scalarField& c,
163  const scalar T,
164  const scalar p,
165  scalar& pf,
166  scalar& cf,
167  label& lRef,
168  scalar& pr,
169  scalar& cr,
170  label& rRef
171 ) const
172 {
173  scalarField c2(nSpecie_, 0.0);
174  for (label i = 0; i < nSpecie_; i++)
175  {
176  c2[i] = max(0.0, c[i]);
177  }
178 
179  const scalar kf = R.kf(p, T, c2);
180  const scalar kr = R.kr(kf, p, T, c2);
181 
182  pf = 1.0;
183  pr = 1.0;
184 
185  const label Nl = R.lhs().size();
186  const label Nr = R.rhs().size();
187 
188  label slRef = 0;
189  lRef = R.lhs()[slRef].index;
190 
191  pf = kf;
192  for (label s = 1; s < Nl; s++)
193  {
194  const label si = R.lhs()[s].index;
195 
196  if (c[si] < c[lRef])
197  {
198  const scalar exp = R.lhs()[slRef].exponent;
199  pf *= pow(max(0.0, c[lRef]), exp);
200  lRef = si;
201  slRef = s;
202  }
203  else
204  {
205  const scalar exp = R.lhs()[s].exponent;
206  pf *= pow(max(0.0, c[si]), exp);
207  }
208  }
209  cf = max(0.0, c[lRef]);
210 
211  {
212  const scalar exp = R.lhs()[slRef].exponent;
213  if (exp < 1.0)
214  {
215  if (cf > SMALL)
216  {
217  pf *= pow(cf, exp - 1.0);
218  }
219  else
220  {
221  pf = 0.0;
222  }
223  }
224  else
225  {
226  pf *= pow(cf, exp - 1.0);
227  }
228  }
229 
230  label srRef = 0;
231  rRef = R.rhs()[srRef].index;
232 
233  // find the matrix element and element position for the rhs
234  pr = kr;
235  for (label s = 1; s < Nr; s++)
236  {
237  const label si = R.rhs()[s].index;
238  if (c[si] < c[rRef])
239  {
240  const scalar exp = R.rhs()[srRef].exponent;
241  pr *= pow(max(0.0, c[rRef]), exp);
242  rRef = si;
243  srRef = s;
244  }
245  else
246  {
247  const scalar exp = R.rhs()[s].exponent;
248  pr *= pow(max(0.0, c[si]), exp);
249  }
250  }
251  cr = max(0.0, c[rRef]);
252 
253  {
254  const scalar exp = R.rhs()[srRef].exponent;
255  if (exp < 1.0)
256  {
257  if (cr>SMALL)
258  {
259  pr *= pow(cr, exp - 1.0);
260  }
261  else
262  {
263  pr = 0.0;
264  }
265  }
266  else
267  {
268  pr *= pow(cr, exp - 1.0);
269  }
270  }
271 
272  return pf*cf - pr*cr;
273 }
274 
275 
276 template<class CompType, class ThermoType>
278 (
279  const scalar time,
280  const scalarField &c,
281  scalarField& dcdt
282 ) const
283 {
284  const scalar T = c[nSpecie_];
285  const scalar p = c[nSpecie_ + 1];
286 
287  dcdt = omega(c, T, p);
288 
289  // constant pressure
290  // dT/dt = ...
291  scalar rho = 0.0;
292  scalar cSum = 0.0;
293  for (label i = 0; i < nSpecie_; i++)
294  {
295  const scalar W = specieThermo_[i].W();
296  cSum += c[i];
297  rho += W*c[i];
298  }
299  scalar cp = 0.0;
300  for (label i=0; i<nSpecie_; i++)
301  {
302  cp += c[i]*specieThermo_[i].cp(p, T);
303  }
304  cp /= rho;
305 
306  scalar dT = 0.0;
307  for (label i = 0; i < nSpecie_; i++)
308  {
309  const scalar hi = specieThermo_[i].ha(p, T);
310  dT += hi*dcdt[i];
311  }
312  dT /= rho*cp;
313 
314  dcdt[nSpecie_] = -dT;
315 
316  // dp/dt = ...
317  dcdt[nSpecie_ + 1] = 0.0;
318 }
319 
320 
321 template<class CompType, class ThermoType>
323 (
324  const scalar t,
325  const scalarField& c,
326  scalarField& dcdt,
327  scalarSquareMatrix& dfdc
328 ) const
329 {
330  const scalar T = c[nSpecie_];
331  const scalar p = c[nSpecie_ + 1];
332 
333  scalarField c2(nSpecie_, 0.0);
334  forAll(c2, i)
335  {
336  c2[i] = max(c[i], 0.0);
337  }
338 
339  for (label i=0; i<nEqns(); i++)
340  {
341  for (label j=0; j<nEqns(); j++)
342  {
343  dfdc[i][j] = 0.0;
344  }
345  }
346 
347  // Length of the first argument must be nSpecie()
348  dcdt = omega(c2, T, p);
349 
350  forAll(reactions_, ri)
351  {
352  const Reaction<ThermoType>& R = reactions_[ri];
353 
354  const scalar kf0 = R.kf(p, T, c2);
355  const scalar kr0 = R.kr(kf0, p, T, c2);
356 
357  forAll(R.lhs(), j)
358  {
359  const label sj = R.lhs()[j].index;
360  scalar kf = kf0;
361  forAll(R.lhs(), i)
362  {
363  const label si = R.lhs()[i].index;
364  const scalar el = R.lhs()[i].exponent;
365  if (i == j)
366  {
367  if (el < 1.0)
368  {
369  if (c2[si] > SMALL)
370  {
371  kf *= el*pow(c2[si] + VSMALL, el - 1.0);
372  }
373  else
374  {
375  kf = 0.0;
376  }
377  }
378  else
379  {
380  kf *= el*pow(c2[si], el - 1.0);
381  }
382  }
383  else
384  {
385  kf *= pow(c2[si], el);
386  }
387  }
388 
389  forAll(R.lhs(), i)
390  {
391  const label si = R.lhs()[i].index;
392  const scalar sl = R.lhs()[i].stoichCoeff;
393  dfdc[si][sj] -= sl*kf;
394  }
395  forAll(R.rhs(), i)
396  {
397  const label si = R.rhs()[i].index;
398  const scalar sr = R.rhs()[i].stoichCoeff;
399  dfdc[si][sj] += sr*kf;
400  }
401  }
402 
403  forAll(R.rhs(), j)
404  {
405  const label sj = R.rhs()[j].index;
406  scalar kr = kr0;
407  forAll(R.rhs(), i)
408  {
409  const label si = R.rhs()[i].index;
410  const scalar er = R.rhs()[i].exponent;
411  if (i == j)
412  {
413  if (er < 1.0)
414  {
415  if (c2[si] > SMALL)
416  {
417  kr *= er*pow(c2[si] + VSMALL, er - 1.0);
418  }
419  else
420  {
421  kr = 0.0;
422  }
423  }
424  else
425  {
426  kr *= er*pow(c2[si], er - 1.0);
427  }
428  }
429  else
430  {
431  kr *= pow(c2[si], er);
432  }
433  }
434 
435  forAll(R.lhs(), i)
436  {
437  const label si = R.lhs()[i].index;
438  const scalar sl = R.lhs()[i].stoichCoeff;
439  dfdc[si][sj] += sl*kr;
440  }
441  forAll(R.rhs(), i)
442  {
443  const label si = R.rhs()[i].index;
444  const scalar sr = R.rhs()[i].stoichCoeff;
445  dfdc[si][sj] -= sr*kr;
446  }
447  }
448  }
449 
450  // Calculate the dcdT elements numerically
451  const scalar delta = 1.0e-3;
452  const scalarField dcdT0(omega(c2, T - delta, p));
453  const scalarField dcdT1(omega(c2, T + delta, p));
454 
455  for (label i = 0; i < nEqns(); i++)
456  {
457  dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
458  }
459 }
460 
461 
462 template<class CompType, class ThermoType>
465 {
466  scalar pf, cf, pr, cr;
467  label lRef, rRef;
468 
469  const volScalarField rho
470  (
471  IOobject
472  (
473  "rho",
474  this->time().timeName(),
475  this->mesh(),
476  IOobject::NO_READ,
477  IOobject::NO_WRITE,
478  false
479  ),
480  this->thermo().rho()
481  );
482 
484  (
485  new volScalarField
486  (
487  IOobject
488  (
489  "tc",
490  this->time().timeName(),
491  this->mesh(),
492  IOobject::NO_READ,
493  IOobject::NO_WRITE,
494  false
495  ),
496  this->mesh(),
497  dimensionedScalar("zero", dimTime, SMALL),
498  zeroGradientFvPatchScalarField::typeName
499  )
500  );
501 
502  scalarField& tc = ttc();
503  const scalarField& T = this->thermo().T();
504  const scalarField& p = this->thermo().p();
505 
506  const label nReaction = reactions_.size();
507 
508  if (this->chemistry_)
509  {
510  forAll(rho, celli)
511  {
512  scalar rhoi = rho[celli];
513  scalar Ti = T[celli];
514  scalar pi = p[celli];
515  scalarField c(nSpecie_);
516  scalar cSum = 0.0;
517 
518  for (label i=0; i<nSpecie_; i++)
519  {
520  scalar Yi = Y_[i][celli];
521  c[i] = rhoi*Yi/specieThermo_[i].W();
522  cSum += c[i];
523  }
524 
525  forAll(reactions_, i)
526  {
527  const Reaction<ThermoType>& R = reactions_[i];
528 
529  omega(R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef);
530 
531  forAll(R.rhs(), s)
532  {
533  scalar sr = R.rhs()[s].stoichCoeff;
534  tc[celli] += sr*pf*cf;
535  }
536  }
537  tc[celli] = nReaction*cSum/tc[celli];
538  }
539  }
540 
541 
542  ttc().correctBoundaryConditions();
543 
544  return ttc;
545 }
546 
547 
548 template<class CompType, class ThermoType>
551 {
553  (
554  new volScalarField
555  (
556  IOobject
557  (
558  "Sh",
559  this->mesh_.time().timeName(),
560  this->mesh_,
561  IOobject::NO_READ,
562  IOobject::NO_WRITE,
563  false
564  ),
565  this->mesh_,
567  zeroGradientFvPatchScalarField::typeName
568  )
569  );
570 
571 
572  if (this->chemistry_)
573  {
574  scalarField& Sh = tSh();
575 
576  forAll(Y_, i)
577  {
578  forAll(Sh, cellI)
579  {
580  const scalar hi = specieThermo_[i].Hc();
581  Sh[cellI] -= hi*RR_[i][cellI];
582  }
583  }
584  }
585 
586  return tSh;
587 }
588 
589 
590 template<class CompType, class ThermoType>
593 {
595  (
596  new volScalarField
597  (
598  IOobject
599  (
600  "dQ",
601  this->mesh_.time().timeName(),
602  this->mesh_,
603  IOobject::NO_READ,
604  IOobject::NO_WRITE,
605  false
606  ),
607  this->mesh_,
608  dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
609  zeroGradientFvPatchScalarField::typeName
610  )
611  );
612 
613  if (this->chemistry_)
614  {
615  volScalarField& dQ = tdQ();
616  dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
617  }
618 
619  return tdQ;
620 }
621 
622 
623 template<class CompType, class ThermoType>
625 {
626  // nEqns = number of species + temperature + pressure
627  return nSpecie_ + 2;
628 }
629 
630 
631 template<class CompType, class ThermoType>
634 (
635  const label reactionI,
636  const label speciei
637 ) const
638 {
639  scalar pf, cf, pr, cr;
640  label lRef, rRef;
641 
642  const volScalarField rho
643  (
644  IOobject
645  (
646  "rho",
647  this->time().timeName(),
648  this->mesh(),
649  IOobject::NO_READ,
650  IOobject::NO_WRITE,
651  false
652  ),
653  this->thermo().rho()
654  );
655 
657  (
659  (
660  IOobject
661  (
662  "RR",
663  this->mesh().time().timeName(),
664  this->mesh(),
665  IOobject::NO_READ,
666  IOobject::NO_WRITE
667  ),
668  this->mesh(),
670  )
671  );
672 
674 
675  const scalarField& T = this->thermo().T();
676  const scalarField& p = this->thermo().p();
677 
678  forAll(rho, celli)
679  {
680  const scalar rhoi = rho[celli];
681  const scalar Ti = T[celli];
682  const scalar pi = p[celli];
683 
684  scalarField c(nSpecie_, 0.0);
685  for (label i=0; i<nSpecie_; i++)
686  {
687  const scalar Yi = Y_[i][celli];
688  c[i] = rhoi*Yi/specieThermo_[i].W();
689  }
690 
691  const scalar w = omegaI
692  (
693  reactionI,
694  c,
695  Ti,
696  pi,
697  pf,
698  cf,
699  lRef,
700  pr,
701  cr,
702  rRef
703  );
704 
705  RR[celli] = w*specieThermo_[speciei].W();
706 
707  }
708 
709  return tRR;
710 }
711 
712 
713 template<class CompType, class ThermoType>
715 {
716  if (!this->chemistry_)
717  {
718  return;
719  }
720 
721  const volScalarField rho
722  (
723  IOobject
724  (
725  "rho",
726  this->time().timeName(),
727  this->mesh(),
728  IOobject::NO_READ,
729  IOobject::NO_WRITE,
730  false
731  ),
732  this->thermo().rho()
733  );
734 
735  const scalarField& T = this->thermo().T();
736  const scalarField& p = this->thermo().p();
737 
738  forAll(rho, celli)
739  {
740  const scalar rhoi = rho[celli];
741  const scalar Ti = T[celli];
742  const scalar pi = p[celli];
743 
744  scalarField c(nSpecie_, 0.0);
745  for (label i=0; i<nSpecie_; i++)
746  {
747  const scalar Yi = Y_[i][celli];
748  c[i] = rhoi*Yi/specieThermo_[i].W();
749  }
750 
751  const scalarField dcdt(omega(c, Ti, pi));
752 
753  for (label i=0; i<nSpecie_; i++)
754  {
755  RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
756  }
757  }
758 }
759 
760 
761 template<class CompType, class ThermoType>
762 template<class DeltaTType>
764 (
765  const DeltaTType& deltaT
766 )
767 {
769 
770  scalar deltaTMin = GREAT;
771 
772  if (!this->chemistry_)
773  {
774  return deltaTMin;
775  }
776 
777  const volScalarField rho
778  (
779  IOobject
780  (
781  "rho",
782  this->time().timeName(),
783  this->mesh(),
784  IOobject::NO_READ,
785  IOobject::NO_WRITE,
786  false
787  ),
788  this->thermo().rho()
789  );
790 
791  const scalarField& T = this->thermo().T();
792  const scalarField& p = this->thermo().p();
793 
794  scalarField c(nSpecie_);
795  scalarField c0(nSpecie_);
796 
797  forAll(rho, celli)
798  {
799  scalar Ti = T[celli];
800 
801  if (Ti > Treact_)
802  {
803  const scalar rhoi = rho[celli];
804  scalar pi = p[celli];
805 
806  for (label i=0; i<nSpecie_; i++)
807  {
808  c[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
809  c0[i] = c[i];
810  }
811 
812  // Initialise time progress
813  scalar timeLeft = deltaT[celli];
814 
815  // Calculate the chemical source terms
816  while (timeLeft > SMALL)
817  {
818  scalar dt = timeLeft;
819  this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
820  timeLeft -= dt;
821  }
822 
823  deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
824 
825  for (label i=0; i<nSpecie_; i++)
826  {
827  RR_[i][celli] =
828  (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
829  }
830  }
831  else
832  {
833  for (label i=0; i<nSpecie_; i++)
834  {
835  RR_[i][celli] = 0;
836  }
837  }
838  }
839 
840  return deltaTMin;
841 }
842 
843 
844 template<class CompType, class ThermoType>
846 (
847  const scalar deltaT
848 )
849 {
850  // Don't allow the time-step to change more than a factor of 2
851  return min
852  (
854  2*deltaT
855  );
856 }
857 
858 
859 template<class CompType, class ThermoType>
861 (
862  const scalarField& deltaT
863 )
864 {
865  return this->solve<scalarField>(deltaT);
866 }
867 
868 
869 template<class CompType, class ThermoType>
871 (
872  scalarField &c,
873  scalar& T,
874  scalar& p,
875  scalar& deltaT,
876  scalar& subDeltaT
877 ) const
878 {
880  (
881  "chemistryModel::solve"
882  "("
883  "scalarField&, "
884  "scalar&, "
885  "scalar&, "
886  "scalar&, "
887  "scalar&"
888  ") const"
889  );
890 }
891 
892 
893 // ************************************************************************* //
virtual label nEqns() const
Number of ODE&#39;s to solve.
virtual tmp< volScalarField > tc() const
Return the chemical time scale.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
scalar Sh
Definition: solveChemistry.H:2
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
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 ))
label nSpecie
rhoEqn solve()
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual scalar kf(const scalar p, const scalar T, const scalarField &c) const
Forward rate constant.
Definition: Reaction.C:456
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil)-psi)*pSat, rhoMin);#1"/home/dm2/henry/OpenFOAM/OpenFOAM-3.0.x/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H"1{alphav=max(min((rho-rholSat)/(rhovSat-rholSat), scalar(1)), scalar(0));alphal=1.0-alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:74
A class for handling words, derived from string.
Definition: word.H:59
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
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:741
DimensionedInternalField & dimensionedInternalField()
Return dimensioned internal field.
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimEnergy
messageStream Info
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
dynamicFvMesh & mesh
const List< specieCoeffs > & rhs() const
Definition: ReactionI.H:59
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
virtual tmp< scalarField > omega(const scalarField &c, const scalar T, const scalar p) const
dc/dt = omega, rate of change in concentration, for each species
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:46
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
basicMultiComponentMixture & composition
Definition: createFields.H:35
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual void calculate()
Calculates the reaction rates.
Foam::reactingMixture.
PtrList< volScalarField > & Y
Definition: createFields.H:36
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual scalar kr(const scalar kfwd, const scalar p, const scalar T, const scalarField &c) const
Reverse rate constant from the given forward rate constant.
Definition: Reaction.C:468
#define forAll(list, i)
Definition: UList.H:421
scalar delta
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
dQ
Definition: YEEqn.H:14
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:53
const List< specieCoeffs > & lhs() const
Definition: ReactionI.H:51
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 and the reference.
const scalar RR
Universal gas constant (default in [J/(kmol K)])
virtual tmp< volScalarField > dQ() const
Return the heat release, i.e. enthalpy/sec [kg/m2/s3].
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
const dimensionedScalar c
Speed of light in a vacuum.
virtual ~chemistryModel()
Destructor.
virtual tmp< DimensionedField< scalar, volMesh > > calculateRR(const label reactionI, const label speciei) const
Return reaction rate of the speciei in reactionI.
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
psiReactionThermo & thermo
Definition: createFields.H:32
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:356
volScalarField scalarField(fieldObject, mesh)
virtual tmp< volScalarField > Sh() const
Return source for enthalpy equation [kg/m/s3].
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
A class representing the concept of a uniform field which stores only the single value and providing ...
Definition: UniformField.H:47
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
word timeName
Definition: getTimeIndex.H:3
const volScalarField & cp