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>
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);
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 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/v2212/OpenFOAM-v2212/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));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...
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)