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>
484 IOobject::NO_REGISTER
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(),
558 IOobject::NO_REGISTER
565 if (this->chemistry_)
573 const scalar hi = specieThermo_[i].Hc();
574 Qdot[celli] -= hi*RR_[i][celli];
577 tQdot.
ref().correctBoundaryConditions();
584 template<
class ReactionThermo,
class ThermoType>
592 scalar pf, cf, pr, cr;
622 const scalar rhoi =
rho[celli];
623 const scalar Ti =
T[celli];
624 const scalar
pi =
p[celli];
626 for (label i=0; i<nSpecie_; i++)
628 const scalar Yi = Y_[i][celli];
629 c_[i] = rhoi*Yi/specieThermo_[i].W();
632 const scalar w = omegaI
646 RR[celli] = w*specieThermo_[si].W();
653 template<
class ReactionThermo,
class ThermoType>
656 if (!this->chemistry_)
669 const scalar rhoi =
rho[celli];
670 const scalar Ti =
T[celli];
671 const scalar
pi =
p[celli];
673 for (label i=0; i<nSpecie_; i++)
675 const scalar Yi = Y_[i][celli];
676 c_[i] = rhoi*Yi/specieThermo_[i].W();
679 omega(c_, Ti,
pi, dcdt_);
681 for (label i=0; i<nSpecie_; i++)
683 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
689 template<
class ReactionThermo,
class ThermoType>
690 template<
class DeltaTType>
693 const DeltaTType& deltaT
698 scalar deltaTMin = GREAT;
700 if (!this->chemistry_)
715 scalar Ti =
T[celli];
719 const scalar rhoi =
rho[celli];
720 scalar
pi =
p[celli];
722 for (label i=0; i<nSpecie_; i++)
724 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
729 scalar timeLeft = deltaT[celli];
732 while (timeLeft > SMALL)
734 scalar dt = timeLeft;
735 this->
solve(c_, Ti,
pi, dt, this->deltaTChem_[celli]);
739 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
741 this->deltaTChem_[celli] =
742 min(this->deltaTChem_[celli], this->deltaTChemMax_);
744 for (label i=0; i<nSpecie_; i++)
747 (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
752 for (label i=0; i<nSpecie_; i++)
763 template<
class ReactionThermo,
class ThermoType>
778 template<
class ReactionThermo,
class ThermoType>
784 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
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
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
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/v2312/OpenFOAM-v2312/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)