33 template<
class ReactionThermo,
class ThermoType>
36 ReactionThermo& thermo
49 (this->
thermo()).speciesData()
53 nReaction_(reactions_.size()),
76 "RR." + Y_[fieldi].
name(),
88 Info<<
"StandardChemistryModel: Number of species = " << nSpecie_
89 <<
" and reactions = " << nReaction_ <<
endl;
95 template<
class ReactionThermo,
class ThermoType>
103 template<
class ReactionThermo,
class ThermoType>
119 R.
omega(p, T, c, dcdt);
124 template<
class ReactionThermo,
class ThermoType>
140 scalar w = R.
omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
145 template<
class ReactionThermo,
class ThermoType>
153 const scalar T = c[nSpecie_];
154 const scalar p = c[nSpecie_ + 1];
158 c_[i] =
max(c[i], 0);
161 omega(c_, T, p, dcdt);
167 for (
label i = 0; i < nSpecie_; i++)
169 const scalar
W = specieThermo_[i].W();
174 for (
label i=0; i<nSpecie_; i++)
176 cp += c_[i]*specieThermo_[i].cp(p, T);
181 for (
label i = 0; i < nSpecie_; i++)
183 const scalar hi = specieThermo_[i].ha(p, T);
188 dcdt[nSpecie_] = -dT;
191 dcdt[nSpecie_ + 1] = 0;
195 template<
class ReactionThermo,
class ThermoType>
204 const scalar T = c[nSpecie_];
205 const scalar p = c[nSpecie_ + 1];
209 c_[i] =
max(c[i], 0);
219 for (
label i = 0; i < nSpecie_; i++)
221 hi[i] = specieThermo_[i].ha(p, T);
222 cpi[i] = specieThermo_[i].cp(p, T);
230 R.
dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd,
false, dummy);
231 R.
dwdT(p, T, c_, omegaI, kfwd, kbwd, J,
false, dummy, nSpecie_);
237 scalar dcpdTMean = 0;
238 for (
label i=0; i<nSpecie_; i++)
240 cpMean += c_[i]*cpi[i];
241 dcpdTMean += c_[i]*specieThermo_[i].dcpdT(p, T);
244 for (
label i=0; i<nSpecie_; i++)
246 dTdt += hi[i]*dcdt[i];
250 for (
label i = 0; i < nSpecie_; i++)
253 for (
label j = 0; j < nSpecie_; j++)
255 J(nSpecie_, i) += hi[j]*J(j, i);
257 J(nSpecie_, i) += cpi[i]*dTdt;
258 J(nSpecie_, i) /= -cpMean;
262 J(nSpecie_, nSpecie_) = 0;
263 for (
label i = 0; i < nSpecie_; i++)
265 J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
267 J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
268 J(nSpecie_, nSpecie_) /= -cpMean;
269 J(nSpecie_, nSpecie_) += dTdt/
T;
273 template<
class ReactionThermo,
class ThermoType>
292 extrapolatedCalculatedFvPatchScalarField::typeName
304 const label nReaction = reactions_.size();
306 scalar pf, cf, pr, cr;
309 if (this->chemistry_)
313 const scalar rhoi = rho[celli];
314 const scalar Ti = T[celli];
315 const scalar
pi = p[celli];
319 for (
label i=0; i<nSpecie_; i++)
321 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
329 R.
omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
333 tc[celli] += R.
rhs()[
s].stoichCoeff*pf*cf;
337 tc[celli] = nReaction*cSum/tc[celli];
341 ttc.
ref().correctBoundaryConditions();
347 template<
class ReactionThermo,
class ThermoType>
358 this->mesh_.time().timeName(),
369 if (this->chemistry_)
377 const scalar hi = specieThermo_[i].Hc();
378 Qdot[celli] -= hi*RR_[i][celli];
387 template<
class ReactionThermo,
class ThermoType>
420 scalar pf, cf, pr, cr;
425 const scalar rhoi = rho[celli];
426 const scalar Ti = T[celli];
427 const scalar
pi = p[celli];
429 for (
label i=0; i<nSpecie_; i++)
431 const scalar Yi = Y_[i][celli];
432 c_[i] = rhoi*Yi/specieThermo_[i].W();
436 const scalar omegai = R.
omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
440 if (si == R.
lhs()[
s].index)
442 RR[celli] -= R.
lhs()[
s].stoichCoeff*omegai;
448 if (si == R.
rhs()[
s].index)
450 RR[celli] += R.
rhs()[
s].stoichCoeff*omegai;
454 RR[celli] *= specieThermo_[si].W();
461 template<
class ReactionThermo,
class ThermoType>
464 if (!this->chemistry_)
477 const scalar rhoi = rho[celli];
478 const scalar Ti = T[celli];
479 const scalar
pi = p[celli];
481 for (
label i=0; i<nSpecie_; i++)
483 const scalar Yi = Y_[i][celli];
484 c_[i] = rhoi*Yi/specieThermo_[i].W();
487 omega(c_, Ti, pi, dcdt_);
489 for (
label i=0; i<nSpecie_; i++)
491 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
497 template<
class ReactionThermo,
class ThermoType>
498 template<
class DeltaTType>
501 const DeltaTType& deltaT
506 scalar deltaTMin = great;
508 if (!this->chemistry_)
523 scalar Ti = T[celli];
527 const scalar rhoi = rho[celli];
528 scalar
pi = p[celli];
530 for (
label i=0; i<nSpecie_; i++)
532 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
537 scalar timeLeft = deltaT[celli];
540 while (timeLeft > small)
542 scalar dt = timeLeft;
543 this->
solve(c_, Ti, pi, dt, this->deltaTChem_[celli]);
547 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
549 this->deltaTChem_[celli] =
550 min(this->deltaTChem_[celli], this->deltaTChemMax_);
552 for (
label i=0; i<nSpecie_; i++)
555 (c_[i] - c0[i])*specieThermo_[i].
W()/deltaT[celli];
560 for (
label i=0; i<nSpecie_; i++)
571 template<
class ReactionThermo,
class ThermoType>
586 template<
class ReactionThermo,
class ThermoType>
592 return this->solve<scalarField>(deltaT);
#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.
Basic chemistry model templated on thermodynamics.
virtual tmp< volScalarField > Qdot() const
Return the heat release rate [kg/m/s3].
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Abstract base class for the systems of ordinary differential equations.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
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.
void dwdT(const scalar p, const scalar T, const scalarField &c, const scalar omegaI, const scalar kfwd, const scalar kbwd, scalarSquareMatrix &J, const bool reduced, const List< label > &c2s, const label indexT) const
Derivative of the net reaction rate for each species involved.
basicSpecieMixture & composition
T & ref() const
Return non-const reference or generate a fatal error.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
rhoReactionThermo & thermo
void omega(const scalar p, const scalar T, const scalarField &c, scalarField &dcdt) const
Net reaction rate for individual species.
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...
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))
virtual void calculate()
Calculates the reaction rates.
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &J) const
Calculate the Jacobian of the system.
void dwdc(const scalar p, const scalar T, const scalarField &c, scalarSquareMatrix &J, scalarField &dcdt, scalar &omegaI, scalar &kfwd, scalar &kbwd, const bool reduced, const List< label > &c2s) const
Derivative of the net reaction rate for each species involved.
virtual ~StandardChemistryModel()
Destructor.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
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-6/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 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.
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const scalarList W(::W(thermo))
PtrList< volScalarField > & Y
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.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
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