StandardChemistryModel.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) 2011-2018 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 "StandardChemistryModel.H"
27 #include "reactingMixture.H"
28 #include "UniformField.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class ReactionThermo, class ThermoType>
35 (
36  ReactionThermo& thermo
37 )
38 :
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_
55  (
56  BasicChemistryModel<ReactionThermo>::template lookupOrDefault<scalar>
57  (
58  "Treact",
59  0
60  )
61  ),
62  RR_(nSpecie_),
63  c_(nSpecie_),
64  dcdt_(nSpecie_)
65 {
66  // Create the fields for the chemistry sources
67  forAll(RR_, fieldi)
68  {
69  RR_.set
70  (
71  fieldi,
73  (
74  IOobject
75  (
76  "RR." + Y_[fieldi].name(),
77  this->mesh().time().timeName(),
78  this->mesh(),
79  IOobject::NO_READ,
80  IOobject::NO_WRITE
81  ),
82  thermo.p().mesh(),
84  )
85  );
86  }
87 
88  Info<< "StandardChemistryModel: Number of species = " << nSpecie_
89  << " and reactions = " << nReaction_ << endl;
90 }
91 
92 
93 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
94 
95 template<class ReactionThermo, class ThermoType>
98 {}
99 
100 
101 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
102 
103 template<class ReactionThermo, class ThermoType>
105 (
106  const scalarField& c,
107  const scalar T,
108  const scalar p,
109  scalarField& dcdt
110 ) const
111 {
112 
113  dcdt = Zero;
114 
115  forAll(reactions_, i)
116  {
117  const Reaction<ThermoType>& R = reactions_[i];
118 
119  R.omega(p, T, c, dcdt);
120  }
121 }
122 
123 
124 template<class ReactionThermo, class ThermoType>
126 (
127  const label index,
128  const scalarField& c,
129  const scalar T,
130  const scalar p,
131  scalar& pf,
132  scalar& cf,
133  label& lRef,
134  scalar& pr,
135  scalar& cr,
136  label& rRef
137 ) const
138 {
139  const Reaction<ThermoType>& R = reactions_[index];
140  scalar w = R.omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
141  return(w);
142 }
143 
144 
145 template<class ReactionThermo, class ThermoType>
147 (
148  const scalar time,
149  const scalarField& c,
150  scalarField& dcdt
151 ) const
152 {
153  const scalar T = c[nSpecie_];
154  const scalar p = c[nSpecie_ + 1];
155 
156  forAll(c_, i)
157  {
158  c_[i] = max(c[i], 0);
159  }
160 
161  omega(c_, T, p, dcdt);
162 
163  // Constant pressure
164  // dT/dt = ...
165  scalar rho = 0;
166  scalar cSum = 0;
167  for (label i = 0; i < nSpecie_; i++)
168  {
169  const scalar W = specieThermo_[i].W();
170  cSum += c_[i];
171  rho += W*c_[i];
172  }
173  scalar cp = 0;
174  for (label i=0; i<nSpecie_; i++)
175  {
176  cp += c_[i]*specieThermo_[i].cp(p, T);
177  }
178  cp /= rho;
179 
180  scalar dT = 0;
181  for (label i = 0; i < nSpecie_; i++)
182  {
183  const scalar hi = specieThermo_[i].ha(p, T);
184  dT += hi*dcdt[i];
185  }
186  dT /= rho*cp;
187 
188  dcdt[nSpecie_] = -dT;
189 
190  // dp/dt = ...
191  dcdt[nSpecie_ + 1] = 0;
192 }
193 
194 
195 template<class ReactionThermo, class ThermoType>
197 (
198  const scalar t,
199  const scalarField& c,
200  scalarField& dcdt,
202 ) const
203 {
204  const scalar T = c[nSpecie_];
205  const scalar p = c[nSpecie_ + 1];
206 
207  forAll(c_, i)
208  {
209  c_[i] = max(c[i], 0);
210  }
211 
212  J = Zero;
213  dcdt = Zero;
214 
215  // To compute the species derivatives of the temperature term,
216  // the enthalpies of the individual species is needed
217  scalarField hi(nSpecie_);
218  scalarField cpi(nSpecie_);
219  for (label i = 0; i < nSpecie_; i++)
220  {
221  hi[i] = specieThermo_[i].ha(p, T);
222  cpi[i] = specieThermo_[i].cp(p, T);
223  }
224  scalar omegaI = 0;
225  List<label> dummy;
226  forAll(reactions_, ri)
227  {
228  const Reaction<ThermoType>& R = reactions_[ri];
229  scalar kfwd, kbwd;
230  R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
231  R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
232  }
233 
234  // The species derivatives of the temperature term are partially computed
235  // while computing dwdc, they are completed hereunder:
236  scalar cpMean = 0;
237  scalar dcpdTMean = 0;
238  for (label i=0; i<nSpecie_; i++)
239  {
240  cpMean += c_[i]*cpi[i]; // J/(m3.K)
241  dcpdTMean += c_[i]*specieThermo_[i].dcpdT(p, T);
242  }
243  scalar dTdt = 0.0;
244  for (label i=0; i<nSpecie_; i++)
245  {
246  dTdt += hi[i]*dcdt[i]; // J/(m3.s)
247  }
248  dTdt /= -cpMean; // K/s
249 
250  for (label i = 0; i < nSpecie_; i++)
251  {
252  J(nSpecie_, i) = 0;
253  for (label j = 0; j < nSpecie_; j++)
254  {
255  J(nSpecie_, i) += hi[j]*J(j, i);
256  }
257  J(nSpecie_, i) += cpi[i]*dTdt; // J/(mol.s)
258  J(nSpecie_, i) /= -cpMean; // K/s / (mol/m3)
259  }
260 
261  // ddT of dTdt
262  J(nSpecie_, nSpecie_) = 0;
263  for (label i = 0; i < nSpecie_; i++)
264  {
265  J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
266  }
267  J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
268  J(nSpecie_, nSpecie_) /= -cpMean;
269  J(nSpecie_, nSpecie_) += dTdt/T;
270 }
271 
272 
273 template<class ReactionThermo, class ThermoType>
276 {
278  (
279  new volScalarField
280  (
281  IOobject
282  (
283  "tc",
284  this->time().timeName(),
285  this->mesh(),
286  IOobject::NO_READ,
287  IOobject::NO_WRITE,
288  false
289  ),
290  this->mesh(),
291  dimensionedScalar("zero", dimTime, small),
292  extrapolatedCalculatedFvPatchScalarField::typeName
293  )
294  );
295 
296  scalarField& tc = ttc.ref();
297 
298  tmp<volScalarField> trho(this->thermo().rho());
299  const scalarField& rho = trho();
300 
301  const scalarField& T = this->thermo().T();
302  const scalarField& p = this->thermo().p();
303 
304  const label nReaction = reactions_.size();
305 
306  scalar pf, cf, pr, cr;
307  label lRef, rRef;
308 
309  if (this->chemistry_)
310  {
311  forAll(rho, celli)
312  {
313  const scalar rhoi = rho[celli];
314  const scalar Ti = T[celli];
315  const scalar pi = p[celli];
316 
317  scalar cSum = 0;
318 
319  for (label i=0; i<nSpecie_; i++)
320  {
321  c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
322  cSum += c_[i];
323  }
324 
325  forAll(reactions_, i)
326  {
327  const Reaction<ThermoType>& R = reactions_[i];
328 
329  R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
330 
331  forAll(R.rhs(), s)
332  {
333  tc[celli] += R.rhs()[s].stoichCoeff*pf*cf;
334  }
335  }
336 
337  tc[celli] = nReaction*cSum/tc[celli];
338  }
339  }
340 
341  ttc.ref().correctBoundaryConditions();
342 
343  return ttc;
344 }
345 
346 
347 template<class ReactionThermo, class ThermoType>
350 {
351  tmp<volScalarField> tQdot
352  (
353  new volScalarField
354  (
355  IOobject
356  (
357  "Qdot",
358  this->mesh_.time().timeName(),
359  this->mesh_,
360  IOobject::NO_READ,
361  IOobject::NO_WRITE,
362  false
363  ),
364  this->mesh_,
366  )
367  );
368 
369  if (this->chemistry_)
370  {
371  scalarField& Qdot = tQdot.ref();
372 
373  forAll(Y_, i)
374  {
375  forAll(Qdot, celli)
376  {
377  const scalar hi = specieThermo_[i].Hc();
378  Qdot[celli] -= hi*RR_[i][celli];
379  }
380  }
381  }
382 
383  return tQdot;
384 }
385 
386 
387 template<class ReactionThermo, class ThermoType>
390 (
391  const label ri,
392  const label si
393 ) const
394 {
396  (
398  (
399  IOobject
400  (
401  "RR",
402  this->mesh().time().timeName(),
403  this->mesh(),
404  IOobject::NO_READ,
405  IOobject::NO_WRITE
406  ),
407  this->mesh(),
409  )
410  );
411 
413 
414  tmp<volScalarField> trho(this->thermo().rho());
415  const scalarField& rho = trho();
416 
417  const scalarField& T = this->thermo().T();
418  const scalarField& p = this->thermo().p();
419 
420  scalar pf, cf, pr, cr;
421  label lRef, rRef;
422 
423  forAll(rho, celli)
424  {
425  const scalar rhoi = rho[celli];
426  const scalar Ti = T[celli];
427  const scalar pi = p[celli];
428 
429  for (label i=0; i<nSpecie_; i++)
430  {
431  const scalar Yi = Y_[i][celli];
432  c_[i] = rhoi*Yi/specieThermo_[i].W();
433  }
434 
435  const Reaction<ThermoType>& R = reactions_[ri];
436  const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
437 
438  forAll(R.lhs(), s)
439  {
440  if (si == R.lhs()[s].index)
441  {
442  RR[celli] -= R.lhs()[s].stoichCoeff*omegai;
443  }
444  }
445 
446  forAll(R.rhs(), s)
447  {
448  if (si == R.rhs()[s].index)
449  {
450  RR[celli] += R.rhs()[s].stoichCoeff*omegai;
451  }
452  }
453 
454  RR[celli] *= specieThermo_[si].W();
455  }
456 
457  return tRR;
458 }
459 
460 
461 template<class ReactionThermo, class ThermoType>
463 {
464  if (!this->chemistry_)
465  {
466  return;
467  }
468 
469  tmp<volScalarField> trho(this->thermo().rho());
470  const scalarField& rho = trho();
471 
472  const scalarField& T = this->thermo().T();
473  const scalarField& p = this->thermo().p();
474 
475  forAll(rho, celli)
476  {
477  const scalar rhoi = rho[celli];
478  const scalar Ti = T[celli];
479  const scalar pi = p[celli];
480 
481  for (label i=0; i<nSpecie_; i++)
482  {
483  const scalar Yi = Y_[i][celli];
484  c_[i] = rhoi*Yi/specieThermo_[i].W();
485  }
486 
487  omega(c_, Ti, pi, dcdt_);
488 
489  for (label i=0; i<nSpecie_; i++)
490  {
491  RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
492  }
493  }
494 }
495 
496 
497 template<class ReactionThermo, class ThermoType>
498 template<class DeltaTType>
500 (
501  const DeltaTType& deltaT
502 )
503 {
505 
506  scalar deltaTMin = great;
507 
508  if (!this->chemistry_)
509  {
510  return deltaTMin;
511  }
512 
513  tmp<volScalarField> trho(this->thermo().rho());
514  const scalarField& rho = trho();
515 
516  const scalarField& T = this->thermo().T();
517  const scalarField& p = this->thermo().p();
518 
519  scalarField c0(nSpecie_);
520 
521  forAll(rho, celli)
522  {
523  scalar Ti = T[celli];
524 
525  if (Ti > Treact_)
526  {
527  const scalar rhoi = rho[celli];
528  scalar pi = p[celli];
529 
530  for (label i=0; i<nSpecie_; i++)
531  {
532  c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
533  c0[i] = c_[i];
534  }
535 
536  // Initialise time progress
537  scalar timeLeft = deltaT[celli];
538 
539  // Calculate the chemical source terms
540  while (timeLeft > small)
541  {
542  scalar dt = timeLeft;
543  this->solve(c_, Ti, pi, dt, this->deltaTChem_[celli]);
544  timeLeft -= dt;
545  }
546 
547  deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
548 
549  this->deltaTChem_[celli] =
550  min(this->deltaTChem_[celli], this->deltaTChemMax_);
551 
552  for (label i=0; i<nSpecie_; i++)
553  {
554  RR_[i][celli] =
555  (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
556  }
557  }
558  else
559  {
560  for (label i=0; i<nSpecie_; i++)
561  {
562  RR_[i][celli] = 0;
563  }
564  }
565  }
566 
567  return deltaTMin;
568 }
569 
570 
571 template<class ReactionThermo, class ThermoType>
573 (
574  const scalar deltaT
575 )
576 {
577  // Don't allow the time-step to change more than a factor of 2
578  return min
579  (
581  2*deltaT
582  );
583 }
584 
585 
586 template<class ReactionThermo, class ThermoType>
588 (
589  const scalarField& deltaT
590 )
591 {
592  return this->solve<scalarField>(deltaT);
593 }
594 
595 
596 // ************************************************************************* //
A class representing the concept of a uniform field which stores only the single value and providing ...
Definition: UniformField.H:47
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
Foam::reactingMixture.
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
Basic chemistry model templated on thermodynamics.
virtual tmp< volScalarField > Qdot() const
Return the heat release rate [kg/m/s3].
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:58
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:46
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 and the reference.
void dwdT(const scalar p, const scalar T, const scalarField &c, const scalar omegaI, const scalar kfwd, const scalar kbwd, scalarSquareMatrix &J, const bool reduced, const List< label > &c2s, const label indexT) const
Derivative of the net reaction rate for each species involved.
Definition: Reaction.C:737
basicSpecieMixture & composition
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
tmp< volScalarField > trho
rhoReactionThermo & thermo
Definition: createFields.H:28
void omega(const scalar p, const scalar T, const scalarField &c, scalarField &dcdt) const
Net reaction rate for individual species.
Definition: Reaction.C:428
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:55
scalar Qdot
Definition: solveChemistry.H:2
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
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))
virtual void calculate()
Calculates the reaction rates.
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &J) const
Calculate the Jacobian of the system.
dynamicFvMesh & mesh
rhoEqn solve()
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:97
void dwdc(const scalar p, const scalar T, const scalarField &c, scalarSquareMatrix &J, scalarField &dcdt, scalar &omegaI, scalar &kfwd, scalar &kbwd, const bool reduced, const List< label > &c2s) const
Derivative of the net reaction rate for each species involved.
Definition: Reaction.C:574
virtual ~StandardChemistryModel()
Destructor.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
const volScalarField & cp
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/ubuntu/OpenFOAM-6/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:72
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const volScalarField & T
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const dimensionSet dimEnergy
virtual tmp< volScalarField::Internal > calculateRR(const label reactionI, const label speciei) const
Return reaction rate of the speciei in reactionI.
virtual tmp< volScalarField > tc() const
Return the chemical time scale.
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const scalarList W(::W(thermo))
PtrList< volScalarField > & Y
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
messageStream Info
const scalar RR
Universal gas constant (default in [J/(kmol K)])
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: ReactionI.H:66
virtual void omega(const scalarField &c, const scalar T, const scalar p, scalarField &dcdt) const
dc/dt = omega, rate of change in concentration, for each species