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>
161 const Reaction<ThermoType>&
R = reactions_[index];
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);
351 const Reaction<ThermoType>&
R = reactions_[ri];
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>
476 IOobject::NO_REGISTER,
489 const label nReaction = reactions_.size();
491 scalar pf, cf, pr, cr;
494 if (this->chemistry_)
498 const scalar rhoi =
rho[celli];
499 const scalar Ti =
T[celli];
500 const scalar
pi =
p[celli];
504 for (label i=0; i<nSpecie_; i++)
506 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
514 omega(
R, c_, Ti,
pi, pf, cf, lRef, pr, cr, rRef);
518 tc[celli] +=
R.rhs()[
s].stoichCoeff*pf*cf;
522 tc[celli] = nReaction*cSum/tc[celli];
526 ttc.ref().correctBoundaryConditions();
532 template<
class ReactionThermo,
class ThermoType>
539 IOobject::NO_REGISTER,
544 if (this->chemistry_)
552 const scalar hi = specieThermo_[i].Hc();
553 Qdot[celli] -= hi*RR_[i][celli];
556 tQdot.ref().correctBoundaryConditions();
563 template<
class ReactionThermo,
class ThermoType>
571 scalar pf, cf, pr, cr;
577 IOobject::NO_REGISTER,
581 auto&
RR = tRR.ref();
591 const scalar rhoi =
rho[celli];
592 const scalar Ti =
T[celli];
593 const scalar
pi =
p[celli];
595 for (label i=0; i<nSpecie_; i++)
597 const scalar Yi = Y_[i][celli];
598 c_[i] = rhoi*Yi/specieThermo_[i].W();
601 const scalar w = omegaI
615 RR[celli] = w*specieThermo_[si].W();
622 template<
class ReactionThermo,
class ThermoType>
625 if (!this->chemistry_)
638 const scalar rhoi =
rho[celli];
639 const scalar Ti =
T[celli];
640 const scalar
pi =
p[celli];
642 for (label i=0; i<nSpecie_; i++)
644 const scalar Yi = Y_[i][celli];
645 c_[i] = rhoi*Yi/specieThermo_[i].W();
648 omega(c_, Ti,
pi, dcdt_);
650 for (label i=0; i<nSpecie_; i++)
652 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
658 template<
class ReactionThermo,
class ThermoType>
659 template<
class DeltaTType>
662 const DeltaTType& deltaT
667 scalar deltaTMin = GREAT;
669 if (!this->chemistry_)
684 scalar Ti =
T[celli];
688 const scalar rhoi =
rho[celli];
689 scalar
pi =
p[celli];
691 for (label i=0; i<nSpecie_; i++)
693 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
698 scalar timeLeft = deltaT[celli];
701 while (timeLeft > SMALL)
703 scalar dt = timeLeft;
704 this->
solve(c_, Ti,
pi, dt, this->deltaTChem_[celli]);
708 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
710 this->deltaTChem_[celli] =
711 min(this->deltaTChem_[celli], this->deltaTChemMax_);
713 for (label i=0; i<nSpecie_; i++)
716 (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
721 for (label i=0; i<nSpecie_; i++)
732 template<
class ReactionThermo,
class ThermoType>
747 template<
class ReactionThermo,
class ThermoType>
753 return this->solve<scalarField>(deltaT);
Basic chemistry model templated on thermodynamics.
virtual tmp< volScalarField > Qdot() const
Return the heat release rate [kg/m/s3].
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.
basicSpecieMixture & composition
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
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...
#define forAll(list, i)
Loop across all elements in list.
psiReactionThermo & thermo
const dimensionSet dimVolume(pow3(dimLength))
virtual void calculate()
Calculates the reaction rates.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
dimensionedScalar exp(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr scalar pi(M_PI)
virtual ~StandardChemistryModel()
Destructor.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
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/chef2/andy/OpenFOAM/release/v2406/OpenFOAM-v2406/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=clamp((rho - rholSat)/(rhovSat - rholSat), zero_one{});alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
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.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
PtrList< volScalarField > & Y
const dimensionedScalar c
Speed of light in a vacuum.
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)
messageStream Info
Information stream (stdout output on master, null elsewhere)
const scalar RR
Universal gas constant: default in [J/(kmol K)].
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
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))
Defines the attributes of an object for which implicit objectRegistry management is supported...
const word extrapolatedCalculatedType
A combined zero-gradient and calculated patch field type.
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
static constexpr const zero Zero
Global zero (0)