33 template<
class CompType,
class ThermoType>
40 CompType(mesh, phaseName),
50 (this->
thermo()).speciesData()
54 nReaction_(reactions_.size()),
55 Treact_(CompType::template lookupOrDefault<scalar>(
"Treact", 0.0)),
70 "RR." + Y_[fieldi].
name(),
82 Info<<
"chemistryModel: Number of species = " << nSpecie_
83 <<
" and reactions = " << nReaction_ <<
endl;
89 template<
class CompType,
class ThermoType>
96 template<
class CompType,
class ThermoType>
105 scalar pf, cf, pr, cr;
114 scalar omegai = omega
116 R, c, T, p, pf, cf, lRef, pr, cr, rRef
122 const scalar sl = R.
lhs()[
s].stoichCoeff;
123 dcdt[si] -= sl*omegai;
129 const scalar sr = R.
rhs()[
s].stoichCoeff;
130 dcdt[si] += sr*omegai;
136 template<
class CompType,
class ThermoType>
152 scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
157 template<
class CompType,
class ThermoType>
172 const scalar kf = R.
kf(p, T, c);
173 const scalar kr = R.
kr(kf, p, T, c);
182 lRef = R.
lhs()[slRef].index;
191 const scalar
exp = R.
lhs()[slRef].exponent;
192 pf *=
pow(
max(0.0, c[lRef]), exp);
198 const scalar
exp = R.
lhs()[
s].exponent;
199 pf *=
pow(
max(0.0, c[si]), exp);
202 cf =
max(0.0, c[lRef]);
205 const scalar
exp = R.
lhs()[slRef].exponent;
210 pf *=
pow(cf, exp - 1.0);
219 pf *=
pow(cf, exp - 1.0);
224 rRef = R.
rhs()[srRef].index;
233 const scalar
exp = R.
rhs()[srRef].exponent;
234 pr *=
pow(
max(0.0, c[rRef]), exp);
240 const scalar
exp = R.
rhs()[
s].exponent;
241 pr *=
pow(
max(0.0, c[si]), exp);
244 cr =
max(0.0, c[rRef]);
247 const scalar
exp = R.
rhs()[srRef].exponent;
252 pr *=
pow(cr, exp - 1.0);
261 pr *=
pow(cr, exp - 1.0);
265 return pf*cf - pr*cr;
269 template<
class CompType,
class ThermoType>
277 const scalar T = c[nSpecie_];
278 const scalar p = c[nSpecie_ + 1];
280 for (
label i = 0; i < nSpecie_; i++)
282 c_[i] =
max(0.0, c[i]);
285 omega(c_, T, p, dcdt);
291 for (
label i = 0; i < nSpecie_; i++)
293 const scalar W = specieThermo_[i].W();
298 for (
label i=0; i<nSpecie_; i++)
300 cp += c_[i]*specieThermo_[i].cp(p, T);
305 for (
label i = 0; i < nSpecie_; i++)
307 const scalar hi = specieThermo_[i].ha(p, T);
312 dcdt[nSpecie_] = -dT;
315 dcdt[nSpecie_ + 1] = 0.0;
319 template<
class CompType,
class ThermoType>
328 const scalar T = c[nSpecie_];
329 const scalar p = c[nSpecie_ + 1];
333 c_[i] =
max(c[i], 0.0);
339 omega(c_, T, p, dcdt);
345 const scalar kf0 = R.
kf(p, T, c_);
346 const scalar kr0 = R.
kr(kf0, p, T, c_);
355 const scalar el = R.
lhs()[i].exponent;
362 kf *= el*
pow(c_[si] + VSMALL, el - 1.0);
371 kf *= el*
pow(c_[si], el - 1.0);
376 kf *=
pow(c_[si], el);
383 const scalar sl = R.
lhs()[i].stoichCoeff;
384 dfdc(si, sj) -= sl*kf;
389 const scalar sr = R.
rhs()[i].stoichCoeff;
390 dfdc(si, sj) += sr*kf;
401 const scalar er = R.
rhs()[i].exponent;
408 kr *= er*
pow(c_[si] + VSMALL, er - 1.0);
417 kr *= er*
pow(c_[si], er - 1.0);
422 kr *=
pow(c_[si], er);
429 const scalar sl = R.
lhs()[i].stoichCoeff;
430 dfdc(si, sj) += sl*kr;
435 const scalar sr = R.
rhs()[i].stoichCoeff;
436 dfdc(si, sj) -= sr*kr;
442 const scalar
delta = 1.0e-3;
444 omega(c_, T + delta, p, dcdt_);
445 for (
label i=0; i<nSpecie_; i++)
447 dfdc(i, nSpecie_) = dcdt_[i];
450 omega(c_, T - delta, p, dcdt_);
451 for (
label i=0; i<nSpecie_; i++)
453 dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/delta;
456 dfdc(nSpecie_, nSpecie_) = 0;
457 dfdc(nSpecie_ + 1, nSpecie_) = 0;
461 template<
class CompType,
class ThermoType>
480 extrapolatedCalculatedFvPatchScalarField::typeName
492 const label nReaction = reactions_.size();
494 scalar pf, cf, pr, cr;
497 if (this->chemistry_)
501 const scalar rhoi = rho[celli];
502 const scalar Ti = T[celli];
503 const scalar
pi = p[celli];
507 for (
label i=0; i<nSpecie_; i++)
509 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
517 omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
521 tc[celli] += R.
rhs()[
s].stoichCoeff*pf*cf;
525 tc[celli] = nReaction*cSum/tc[celli];
529 ttc.
ref().correctBoundaryConditions();
535 template<
class CompType,
class ThermoType>
546 this->mesh_.time().timeName(),
557 if (this->chemistry_)
565 const scalar hi = specieThermo_[i].Hc();
566 Qdot[celli] -= hi*RR_[i][celli];
575 template<
class CompType,
class ThermoType>
583 scalar pf, cf, pr, cr;
613 const scalar rhoi = rho[celli];
614 const scalar Ti = T[celli];
615 const scalar
pi = p[celli];
617 for (
label i=0; i<nSpecie_; i++)
619 const scalar Yi = Y_[i][celli];
620 c_[i] = rhoi*Yi/specieThermo_[i].W();
623 const scalar w = omegaI
637 RR[celli] = w*specieThermo_[si].W();
644 template<
class CompType,
class ThermoType>
647 if (!this->chemistry_)
660 const scalar rhoi = rho[celli];
661 const scalar Ti = T[celli];
662 const scalar
pi = p[celli];
664 for (
label i=0; i<nSpecie_; i++)
666 const scalar Yi = Y_[i][celli];
667 c_[i] = rhoi*Yi/specieThermo_[i].W();
670 omega(c_, Ti, pi, dcdt_);
672 for (
label i=0; i<nSpecie_; i++)
674 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
680 template<
class CompType,
class ThermoType>
681 template<
class DeltaTType>
684 const DeltaTType& deltaT
689 scalar deltaTMin = GREAT;
691 if (!this->chemistry_)
706 scalar Ti = T[celli];
710 const scalar rhoi = rho[celli];
711 scalar
pi = p[celli];
713 for (
label i=0; i<nSpecie_; i++)
715 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
720 scalar timeLeft = deltaT[celli];
723 while (timeLeft > SMALL)
725 scalar dt = timeLeft;
726 this->
solve(c_, Ti, pi, dt, this->deltaTChem_[celli]);
730 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
732 for (
label i=0; i<nSpecie_; i++)
735 (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
740 for (
label i=0; i<nSpecie_; i++)
751 template<
class CompType,
class ThermoType>
766 template<
class CompType,
class ThermoType>
772 return this->solve<scalarField>(deltaT);
virtual tmp< volScalarField::Internal > calculateRR(const label reactionI, const label speciei) const
Return reaction rate of the speciei in reactionI.
virtual scalar kr(const scalar kfwd, const scalar p, const scalar T, const scalarField &c) const
Reverse rate constant from the given forward rate constant.
virtual scalar kf(const scalar p, const scalar T, const scalarField &c) const
Forward rate constant.
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.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const List< specieCoeffs > & lhs() const
Abstract base class for the systems of ordinary differential equations.
virtual tmp< volScalarField > tc() const
Return the chemical time scale.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
basicMultiComponentMixture & composition
T & ref() const
Return non-const reference or generate a fatal error.
virtual ~chemistryModel()
Destructor.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
const Time & time() const
Return the top-level database.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
psiReactionThermo & thermo
const dimensionSet dimVolume(pow3(dimLength))
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
dimensionedScalar exp(const dimensionedScalar &ds)
A class for handling words, derived from string.
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/ubuntu/OpenFOAM-5.0/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()
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
const dimensionSet dimEnergy
virtual void calculate()
Calculates the reaction rates.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
PtrList< volScalarField > & Y
Mesh data needed to do the Finite Volume discretisation.
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)
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
const scalar RR
Universal gas constant (default in [J/(kmol K)])
virtual tmp< volScalarField > Qdot() const
Return the heat release rate [kg/m/s3].
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
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
const List< specieCoeffs > & rhs() const