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);
300 for (label i = 0; i < nSpecie_; i++)
302 const scalar W = specieThermo_[i].W();
306 for (label i=0; i<nSpecie_; i++)
308 cp += c_[i]*specieThermo_[i].cp(
p,
T);
313 for (label i = 0; i < nSpecie_; i++)
315 const scalar hi = specieThermo_[i].ha(
p,
T);
320 dcdt[nSpecie_] = -dT;
323 dcdt[nSpecie_ + 1] = 0.0;
327 template<
class ReactionThermo,
class ThermoType>
336 const scalar
T =
c[nSpecie_];
337 const scalar
p =
c[nSpecie_ + 1];
341 c_[i] =
max(
c[i], 0.0);
347 omega(c_,
T,
p, dcdt);
353 const scalar kf0 =
R.kf(
p,
T, c_);
354 const scalar kr0 =
R.kr(kf0,
p,
T, c_);
358 const label sj =
R.lhs()[j].index;
362 const label si =
R.lhs()[i].index;
363 const scalar el =
R.lhs()[i].exponent;
370 kf *= el*
pow(c_[si], el - 1.0);
379 kf *= el*
pow(c_[si], el - 1.0);
384 kf *=
pow(c_[si], el);
390 const label si =
R.lhs()[i].index;
391 const scalar sl =
R.lhs()[i].stoichCoeff;
392 dfdc(si, sj) -= sl*kf;
396 const label si =
R.rhs()[i].index;
397 const scalar sr =
R.rhs()[i].stoichCoeff;
398 dfdc(si, sj) += sr*kf;
404 const label sj =
R.rhs()[j].index;
408 const label si =
R.rhs()[i].index;
409 const scalar er =
R.rhs()[i].exponent;
416 kr *= er*
pow(c_[si], er - 1.0);
425 kr *= er*
pow(c_[si], er - 1.0);
430 kr *=
pow(c_[si], er);
436 const label si =
R.lhs()[i].index;
437 const scalar sl =
R.lhs()[i].stoichCoeff;
438 dfdc(si, sj) += sl*kr;
442 const label si =
R.rhs()[i].index;
443 const scalar sr =
R.rhs()[i].stoichCoeff;
444 dfdc(si, sj) -= sr*kr;
450 const scalar
delta = 1.0e-3;
452 omega(c_,
T +
delta,
p, dcdt_);
453 for (label i=0; i<nSpecie_; i++)
455 dfdc(i, nSpecie_) = dcdt_[i];
458 omega(c_,
T -
delta,
p, dcdt_);
459 for (label i=0; i<nSpecie_; i++)
461 dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/
delta;
464 dfdc(nSpecie_, nSpecie_) = 0;
465 dfdc(nSpecie_ + 1, nSpecie_) = 0;
469 template<
class ReactionThermo,
class ThermoType>
488 extrapolatedCalculatedFvPatchScalarField::typeName
500 const label nReaction = reactions_.size();
502 scalar pf, cf, pr, cr;
505 if (this->chemistry_)
509 const scalar rhoi =
rho[celli];
510 const scalar Ti =
T[celli];
511 const scalar
pi =
p[celli];
515 for (label i=0; i<nSpecie_; i++)
517 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
525 omega(
R, c_, Ti,
pi, pf, cf, lRef, pr, cr, rRef);
529 tc[celli] +=
R.rhs()[
s].stoichCoeff*pf*cf;
533 tc[celli] = nReaction*cSum/tc[celli];
537 ttc.
ref().correctBoundaryConditions();
543 template<
class ReactionThermo,
class ThermoType>
554 this->mesh_.time().timeName(),
565 if (this->chemistry_)
573 const scalar hi = specieThermo_[i].Hc();
574 Qdot[celli] -= hi*RR_[i][celli];
583 template<
class ReactionThermo,
class ThermoType>
591 scalar pf, cf, pr, cr;
621 const scalar rhoi =
rho[celli];
622 const scalar Ti =
T[celli];
623 const scalar
pi =
p[celli];
625 for (label i=0; i<nSpecie_; i++)
627 const scalar Yi = Y_[i][celli];
628 c_[i] = rhoi*Yi/specieThermo_[i].W();
631 const scalar w = omegaI
645 RR[celli] = w*specieThermo_[si].W();
652 template<
class ReactionThermo,
class ThermoType>
655 if (!this->chemistry_)
668 const scalar rhoi =
rho[celli];
669 const scalar Ti =
T[celli];
670 const scalar
pi =
p[celli];
672 for (label i=0; i<nSpecie_; i++)
674 const scalar Yi = Y_[i][celli];
675 c_[i] = rhoi*Yi/specieThermo_[i].W();
678 omega(c_, Ti,
pi, dcdt_);
680 for (label i=0; i<nSpecie_; i++)
682 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
688 template<
class ReactionThermo,
class ThermoType>
689 template<
class DeltaTType>
692 const DeltaTType& deltaT
697 scalar deltaTMin = GREAT;
699 if (!this->chemistry_)
714 scalar Ti =
T[celli];
718 const scalar rhoi =
rho[celli];
719 scalar
pi =
p[celli];
721 for (label i=0; i<nSpecie_; i++)
723 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
728 scalar timeLeft = deltaT[celli];
731 while (timeLeft > SMALL)
733 scalar dt = timeLeft;
734 this->
solve(c_, Ti,
pi, dt, this->deltaTChem_[celli]);
738 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
740 this->deltaTChem_[celli] =
741 min(this->deltaTChem_[celli], this->deltaTChemMax_);
743 for (label i=0; i<nSpecie_; i++)
746 (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
751 for (label i=0; i<nSpecie_; i++)
762 template<
class ReactionThermo,
class ThermoType>
777 template<
class ReactionThermo,
class ThermoType>
783 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 (stdout output on master, null elsewhere)
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.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
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...