36template<
class ReactionThermo,
class ThermoType>
52 (this->
thermo()).speciesData()
56 nReaction_(reactions_.size()),
91 Info<<
"StandardChemistryModel: Number of species = " <<
nSpecie_
98template<
class ReactionThermo,
class ThermoType>
106template<
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;
146template<
class ReactionThermo,
class ThermoType>
162 scalar w = omega(
R, c,
T,
p, pf, cf, lRef, pr, cr, rRef);
167template<
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;
279template<
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;
327template<
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;
469template<
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();
543template<
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];
583template<
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();
652template<
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();
688template<
class ReactionThermo,
class ThermoType>
689template<
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++)
762template<
class ReactionThermo,
class ThermoType>
777template<
class ReactionThermo,
class ThermoType>
783 return this->solve<scalarField>(deltaT);
#define R(A, B, C, D, E, F, K, M)
Basic chemistry model templated on thermodynamics.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const Time & time() const
Return Time associated with the objectRegistry.
Abstract base class for the systems of ordinary differential equations.
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
virtual tmp< volScalarField > tc() const
Return the chemical time scale.
virtual tmp< volScalarField::Internal > calculateRR(const label reactionI, const label speciei) const
Return reaction rate of the speciei in reactionI.
label nSpecie_
Number of species.
PtrList< volScalarField > & Y_
Reference to the field of specie mass fractions.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
virtual ~StandardChemistryModel()
Destructor.
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.
PtrList< volScalarField::Internal > RR_
List of reaction rate per specie [kg/m3/s].
label nReaction_
Number of reactions.
virtual tmp< volScalarField > Qdot() const
Return the heat release rate [kg/m/s3].
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
virtual void calculate()
Calculates the reaction rates.
const word & name() const
const fvMesh & mesh() const
Return const access to the mesh database.
const volScalarField & omega() const
Access functions.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
A class for managing temporary objects.
basicSpecieMixture & composition
PtrList< volScalarField > & Y
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))
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimEnergy
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVolume(pow3(dimLength))
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
tmp< volScalarField > trho
const volScalarField & cp
#define forAll(list, i)
Loop across all elements in list.