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)),
68 "RR." + Y_[fieldi].
name(),
80 Info<<
"chemistryModel: Number of species = " << nSpecie_
81 <<
" and reactions = " << nReaction_ <<
endl;
87 template<
class CompType,
class ThermoType>
94 template<
class CompType,
class ThermoType>
103 scalar pf, cf, pr, cr;
113 scalar omegai = omega
115 R, c, T, p, pf, cf, lRef, pr, cr, rRef
121 const scalar sl = R.
lhs()[
s].stoichCoeff;
128 const scalar sr = R.
rhs()[
s].stoichCoeff;
137 template<
class CompType,
class ThermoType>
154 scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
159 template<
class CompType,
class ThermoType>
175 for (
label i = 0; i < nSpecie_; i++)
177 c2[i] =
max(0.0, c[i]);
180 const scalar kf = R.
kf(p, T, c2);
181 const scalar kr = R.
kr(kf, p, T, c2);
190 lRef = R.
lhs()[slRef].index;
199 const scalar
exp = R.
lhs()[slRef].exponent;
200 pf *=
pow(
max(0.0, c[lRef]), exp);
206 const scalar
exp = R.
lhs()[
s].exponent;
207 pf *=
pow(
max(0.0, c[si]), exp);
210 cf =
max(0.0, c[lRef]);
213 const scalar
exp = R.
lhs()[slRef].exponent;
218 pf *=
pow(cf, exp - 1.0);
227 pf *=
pow(cf, exp - 1.0);
232 rRef = R.
rhs()[srRef].index;
241 const scalar
exp = R.
rhs()[srRef].exponent;
242 pr *=
pow(
max(0.0, c[rRef]), exp);
248 const scalar
exp = R.
rhs()[
s].exponent;
249 pr *=
pow(
max(0.0, c[si]), exp);
252 cr =
max(0.0, c[rRef]);
255 const scalar
exp = R.
rhs()[srRef].exponent;
260 pr *=
pow(cr, exp - 1.0);
269 pr *=
pow(cr, exp - 1.0);
273 return pf*cf - pr*cr;
277 template<
class CompType,
class ThermoType>
285 const scalar T = c[nSpecie_];
286 const scalar p = c[nSpecie_ + 1];
288 dcdt = omega(c, T, p);
294 for (
label i = 0; i < nSpecie_; i++)
296 const scalar W = specieThermo_[i].W();
301 for (
label i=0; i<nSpecie_; i++)
303 cp += c[i]*specieThermo_[i].cp(p, T);
308 for (
label i = 0; i < nSpecie_; i++)
310 const scalar hi = specieThermo_[i].ha(p, T);
315 dcdt[nSpecie_] = -dT;
318 dcdt[nSpecie_ + 1] = 0.0;
322 template<
class CompType,
class ThermoType>
331 const scalar T = c[nSpecie_];
332 const scalar p = c[nSpecie_ + 1];
337 c2[i] =
max(c[i], 0.0);
340 for (
label i=0; i<nEqns(); i++)
342 for (
label j=0; j<nEqns(); j++)
349 dcdt = omega(c2, T, p);
355 const scalar kf0 = R.
kf(p, T, c2);
356 const scalar kr0 = R.
kr(kf0, p, T, c2);
365 const scalar el = R.
lhs()[i].exponent;
372 kf *= el*
pow(c2[si] + VSMALL, el - 1.0);
381 kf *= el*
pow(c2[si], el - 1.0);
386 kf *=
pow(c2[si], el);
393 const scalar sl = R.
lhs()[i].stoichCoeff;
394 dfdc[si][sj] -= sl*kf;
399 const scalar sr = R.
rhs()[i].stoichCoeff;
400 dfdc[si][sj] += sr*kf;
411 const scalar er = R.
rhs()[i].exponent;
418 kr *= er*
pow(c2[si] + VSMALL, er - 1.0);
427 kr *= er*
pow(c2[si], er - 1.0);
432 kr *=
pow(c2[si], er);
439 const scalar sl = R.
lhs()[i].stoichCoeff;
440 dfdc[si][sj] += sl*kr;
445 const scalar sr = R.
rhs()[i].stoichCoeff;
446 dfdc[si][sj] -= sr*kr;
452 const scalar
delta = 1.0e-3;
456 for (
label i = 0; i < nEqns(); i++)
458 dfdc[i][
nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
463 template<
class CompType,
class ThermoType>
467 scalar pf, cf, pr, cr;
499 extrapolatedCalculatedFvPatchScalarField::typeName
507 const label nReaction = reactions_.size();
509 if (this->chemistry_)
513 scalar rhoi = rho[celli];
514 scalar Ti = T[celli];
515 scalar
pi = p[celli];
519 for (
label i=0; i<nSpecie_; i++)
521 scalar Yi = Y_[i][celli];
522 c[i] = rhoi*Yi/specieThermo_[i].W();
530 omega(R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef);
534 scalar sr = R.
rhs()[
s].stoichCoeff;
535 tc[celli] += sr*pf*cf;
538 tc[celli] = nReaction*cSum/tc[celli];
543 ttc.
ref().correctBoundaryConditions();
549 template<
class CompType,
class ThermoType>
560 this->mesh_.time().timeName(),
571 if (this->chemistry_)
579 const scalar hi = specieThermo_[i].Hc();
580 Sh[celli] -= hi*RR_[i][celli];
589 template<
class CompType,
class ThermoType>
600 this->mesh_.time().timeName(),
611 if (this->chemistry_)
614 dQ.
ref() = this->mesh_.V()*
Sh()();
621 template<
class CompType,
class ThermoType>
629 template<
class CompType,
class ThermoType>
633 const label reactionI,
637 scalar pf, cf, pr, cr;
678 const scalar rhoi = rho[celli];
679 const scalar Ti = T[celli];
680 const scalar
pi = p[celli];
683 for (
label i=0; i<nSpecie_; i++)
685 const scalar Yi = Y_[i][celli];
686 c[i] = rhoi*Yi/specieThermo_[i].W();
689 const scalar w = omegaI
703 RR[celli] = w*specieThermo_[speciei].W();
711 template<
class CompType,
class ThermoType>
714 if (!this->chemistry_)
738 const scalar rhoi = rho[celli];
739 const scalar Ti = T[celli];
740 const scalar
pi = p[celli];
743 for (
label i=0; i<nSpecie_; i++)
745 const scalar Yi = Y_[i][celli];
746 c[i] = rhoi*Yi/specieThermo_[i].W();
751 for (
label i=0; i<nSpecie_; i++)
753 RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
759 template<
class CompType,
class ThermoType>
760 template<
class DeltaTType>
763 const DeltaTType& deltaT
768 scalar deltaTMin = GREAT;
770 if (!this->chemistry_)
797 scalar Ti = T[celli];
801 const scalar rhoi = rho[celli];
802 scalar
pi = p[celli];
804 for (
label i=0; i<nSpecie_; i++)
806 c[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
811 scalar timeLeft = deltaT[celli];
814 while (timeLeft > SMALL)
816 scalar dt = timeLeft;
817 this->
solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
821 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
823 for (
label i=0; i<nSpecie_; i++)
826 (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
831 for (
label i=0; i<nSpecie_; i++)
842 template<
class CompType,
class ThermoType>
857 template<
class CompType,
class ThermoType>
863 return this->solve<scalarField>(deltaT);
867 template<
class CompType,
class ThermoType>
const List< specieCoeffs > & lhs() const
virtual tmp< scalarField > omega(const scalarField &c, const scalar T, const scalar p) const
dc/dt = omega, rate of change in concentration, for each species
#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.
Abstract base class for the systems of ordinary differential equations.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
basicMultiComponentMixture & composition
virtual tmp< volScalarField > Sh() const
Return source for enthalpy equation [kg/m/s3].
virtual ~chemistryModel()
Destructor.
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Ostream & endl(Ostream &os)
Add newline and flush stream.
const List< specieCoeffs > & rhs() const
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
virtual scalar kf(const scalar p, const scalar T, const scalarField &c) const
Forward rate constant.
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.
virtual tmp< volScalarField > dQ() const
Return the heat release, i.e. enthalpy/sec [kg/m2/s3].
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
virtual tmp< DimensionedField< scalar, volMesh > > 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.
volScalarField scalarField(fieldObject, mesh)
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-4.1/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()
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.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
const dimensionSet dimEnergy
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.
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
Internal & ref()
Return a reference to the dimensioned internal field.
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.
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)
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.
T & ref() const
Return non-const reference or generate a fatal error.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
virtual label nEqns() const
Number of ODE's to solve.
const Time & time() const
Return the top-level database.