Go to the documentation of this file.
36 template<
class ReactionThermo,
class ThermoType>
52 (this->
thermo()).speciesData()
56 nReaction_(reactions_.size()),
79 "RR." + Y_[fieldi].
name(),
91 Info<<
"StandardChemistryModel: Number of species = " << nSpecie_
92 <<
" and reactions = " << nReaction_ <<
endl;
98 template<
class ReactionThermo,
class ThermoType>
106 template<
class ReactionThermo,
class ThermoType>
115 scalar pf, cf, pr, cr;
124 scalar omegai = omega
126 R,
c,
T,
p, pf, cf, lRef, pr, cr, rRef
131 const label si =
R.lhs()[
s].index;
132 const scalar sl =
R.lhs()[
s].stoichCoeff;
133 dcdt[si] -= sl*omegai;
138 const label si =
R.rhs()[
s].index;
139 const scalar sr =
R.rhs()[
s].stoichCoeff;
140 dcdt[si] += sr*omegai;
146 template<
class ReactionThermo,
class ThermoType>
162 scalar w = omega(
R,
c,
T,
p, pf, cf, lRef, pr, cr, rRef);
167 template<
class ReactionThermo,
class ThermoType>
182 const scalar kf =
R.kf(
p,
T,
c);
183 const scalar kr =
R.kr(kf,
p,
T,
c);
188 const label Nl =
R.lhs().size();
189 const label Nr =
R.rhs().size();
192 lRef =
R.lhs()[slRef].index;
195 for (label
s = 1;
s < Nl;
s++)
197 const label si =
R.lhs()[
s].index;
201 const scalar
exp =
R.lhs()[slRef].exponent;
208 const scalar
exp =
R.lhs()[
s].exponent;
212 cf =
max(
c[lRef], 0.0);
215 const scalar
exp =
R.lhs()[slRef].exponent;
234 rRef =
R.rhs()[srRef].index;
238 for (label
s = 1;
s < Nr;
s++)
240 const label si =
R.rhs()[
s].index;
243 const scalar
exp =
R.rhs()[srRef].exponent;
250 const scalar
exp =
R.rhs()[
s].exponent;
254 cr =
max(
c[rRef], 0.0);
257 const scalar
exp =
R.rhs()[srRef].exponent;
275 return pf*cf - pr*cr;
279 template<
class ReactionThermo,
class ThermoType>
287 const scalar
T =
c[nSpecie_];
288 const scalar
p =
c[nSpecie_ + 1];
292 c_[i] =
max(
c[i], 0.0);
295 omega(c_,
T,
p, dcdt);
301 for (label i = 0; i < nSpecie_; i++)
303 const scalar W = specieThermo_[i].W();
308 for (label i=0; i<nSpecie_; i++)
310 cp += c_[i]*specieThermo_[i].cp(
p,
T);
315 for (label i = 0; i < nSpecie_; i++)
317 const scalar hi = specieThermo_[i].ha(
p,
T);
322 dcdt[nSpecie_] = -dT;
325 dcdt[nSpecie_ + 1] = 0.0;
329 template<
class ReactionThermo,
class ThermoType>
338 const scalar
T =
c[nSpecie_];
339 const scalar
p =
c[nSpecie_ + 1];
343 c_[i] =
max(
c[i], 0.0);
349 omega(c_,
T,
p, dcdt);
355 const scalar kf0 =
R.kf(
p,
T, c_);
356 const scalar kr0 =
R.kr(kf0,
p,
T, c_);
360 const label sj =
R.lhs()[j].index;
364 const label si =
R.lhs()[i].index;
365 const scalar el =
R.lhs()[i].exponent;
372 kf *= el*
pow(c_[si], el - 1.0);
381 kf *= el*
pow(c_[si], el - 1.0);
386 kf *=
pow(c_[si], el);
392 const label si =
R.lhs()[i].index;
393 const scalar sl =
R.lhs()[i].stoichCoeff;
394 dfdc(si, sj) -= sl*kf;
398 const label si =
R.rhs()[i].index;
399 const scalar sr =
R.rhs()[i].stoichCoeff;
400 dfdc(si, sj) += sr*kf;
406 const label sj =
R.rhs()[j].index;
410 const label si =
R.rhs()[i].index;
411 const scalar er =
R.rhs()[i].exponent;
418 kr *= er*
pow(c_[si], er - 1.0);
427 kr *= er*
pow(c_[si], er - 1.0);
432 kr *=
pow(c_[si], er);
438 const label si =
R.lhs()[i].index;
439 const scalar sl =
R.lhs()[i].stoichCoeff;
440 dfdc(si, sj) += sl*kr;
444 const label si =
R.rhs()[i].index;
445 const scalar sr =
R.rhs()[i].stoichCoeff;
446 dfdc(si, sj) -= sr*kr;
452 const scalar
delta = 1.0e-3;
454 omega(c_,
T +
delta,
p, dcdt_);
455 for (label i=0; i<nSpecie_; i++)
457 dfdc(i, nSpecie_) = dcdt_[i];
460 omega(c_,
T -
delta,
p, dcdt_);
461 for (label i=0; i<nSpecie_; i++)
463 dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/
delta;
466 dfdc(nSpecie_, nSpecie_) = 0;
467 dfdc(nSpecie_ + 1, nSpecie_) = 0;
471 template<
class ReactionThermo,
class ThermoType>
490 extrapolatedCalculatedFvPatchScalarField::typeName
502 const label nReaction = reactions_.size();
504 scalar pf, cf, pr, cr;
507 if (this->chemistry_)
511 const scalar rhoi =
rho[celli];
512 const scalar Ti =
T[celli];
513 const scalar
pi =
p[celli];
517 for (label i=0; i<nSpecie_; i++)
519 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
527 omega(
R, c_, Ti,
pi, pf, cf, lRef, pr, cr, rRef);
531 tc[celli] +=
R.rhs()[
s].stoichCoeff*pf*cf;
535 tc[celli] = nReaction*cSum/tc[celli];
539 ttc.
ref().correctBoundaryConditions();
545 template<
class ReactionThermo,
class ThermoType>
556 this->mesh_.time().timeName(),
567 if (this->chemistry_)
575 const scalar hi = specieThermo_[i].Hc();
576 Qdot[celli] -= hi*RR_[i][celli];
585 template<
class ReactionThermo,
class ThermoType>
593 scalar pf, cf, pr, cr;
623 const scalar rhoi =
rho[celli];
624 const scalar Ti =
T[celli];
625 const scalar
pi =
p[celli];
627 for (label i=0; i<nSpecie_; i++)
629 const scalar Yi = Y_[i][celli];
630 c_[i] = rhoi*Yi/specieThermo_[i].W();
633 const scalar w = omegaI
647 RR[celli] = w*specieThermo_[si].W();
654 template<
class ReactionThermo,
class ThermoType>
657 if (!this->chemistry_)
670 const scalar rhoi =
rho[celli];
671 const scalar Ti =
T[celli];
672 const scalar
pi =
p[celli];
674 for (label i=0; i<nSpecie_; i++)
676 const scalar Yi = Y_[i][celli];
677 c_[i] = rhoi*Yi/specieThermo_[i].W();
680 omega(c_, Ti,
pi, dcdt_);
682 for (label i=0; i<nSpecie_; i++)
684 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
690 template<
class ReactionThermo,
class ThermoType>
691 template<
class DeltaTType>
694 const DeltaTType& deltaT
699 scalar deltaTMin = GREAT;
701 if (!this->chemistry_)
716 scalar Ti =
T[celli];
720 const scalar rhoi =
rho[celli];
721 scalar
pi =
p[celli];
723 for (label i=0; i<nSpecie_; i++)
725 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
730 scalar timeLeft = deltaT[celli];
733 while (timeLeft > SMALL)
735 scalar dt = timeLeft;
736 this->
solve(c_, Ti,
pi, dt, this->deltaTChem_[celli]);
740 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
742 this->deltaTChem_[celli] =
743 min(this->deltaTChem_[celli], this->deltaTChemMax_);
745 for (label i=0; i<nSpecie_; i++)
748 (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
753 for (label i=0; i<nSpecie_; i++)
764 template<
class ReactionThermo,
class ThermoType>
779 template<
class ReactionThermo,
class ThermoType>
785 return this->solve<scalarField>(deltaT);
const scalar RR
Universal gas constant: default in [J/(kmol K)].
Defines the attributes of an object for which implicit objectRegistry management is supported,...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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 ~StandardChemistryModel()
Destructor.
const volScalarField & cp
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimEnergy
psiReactionThermo & thermo
Basic chemistry model templated on thermodynamics.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar exp(const dimensionedScalar &ds)
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
basicSpecieMixture & composition
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
#define forAll(list, i)
Loop across all elements in list.
tmp< volScalarField > trho
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
#define R(A, B, C, D, E, F, K, M)
messageStream Info
Information stream (uses stdout - output is on the master only)
word name(const complex &c)
Return string representation of complex.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual void calculate()
Calculates the reaction rates.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
virtual tmp< volScalarField > tc() const
Return the chemical time scale.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
PtrList< volScalarField > & Y
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
constexpr scalar pi(M_PI)
Abstract base class for the systems of ordinary differential equations.
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 dimensionedScalar c
Speed of light in a vacuum.
virtual tmp< volScalarField > Qdot() const
Return the heat release rate [kg/m/s3].
const dimensionSet dimVolume(pow3(dimLength))
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
virtual tmp< volScalarField::Internal > calculateRR(const label reactionI, const label speciei) const
Return reaction rate of the speciei in reactionI.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...