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>
284 extrapolatedCalculatedFvPatchScalarField::typeName
296 const label nReaction = reactions_.size();
298 scalar pf, cf, pr, cr;
301 if (this->chemistry_)
305 const scalar rhoi = rho[celli];
306 const scalar Ti = T[celli];
307 const scalar
pi = p[celli];
311 for (
label i=0; i<nSpecie_; i++)
313 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
321 R.
omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
325 tc[celli] += R.
rhs()[
s].stoichCoeff*pf*cf;
329 tc[celli] = nReaction*cSum/tc[celli];
333 ttc.
ref().correctBoundaryConditions();
339 template<
class ReactionThermo,
class ThermoType>
353 if (this->chemistry_)
361 const scalar hi = specieThermo_[i].Hc();
362 Qdot[celli] -= hi*RR_[i][celli];
371 template<
class ReactionThermo,
class ThermoType>
397 scalar pf, cf, pr, cr;
402 const scalar rhoi = rho[celli];
403 const scalar Ti = T[celli];
404 const scalar
pi = p[celli];
406 for (
label i=0; i<nSpecie_; i++)
408 const scalar Yi = Y_[i][celli];
409 c_[i] = rhoi*Yi/specieThermo_[i].W();
413 const scalar omegai = R.
omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
417 if (si == R.
lhs()[
s].index)
419 RR[celli] -= R.
lhs()[
s].stoichCoeff*omegai;
425 if (si == R.
rhs()[
s].index)
427 RR[celli] += R.
rhs()[
s].stoichCoeff*omegai;
431 RR[celli] *= specieThermo_[si].W();
438 template<
class ReactionThermo,
class ThermoType>
441 if (!this->chemistry_)
454 const scalar rhoi = rho[celli];
455 const scalar Ti = T[celli];
456 const scalar
pi = p[celli];
458 for (
label i=0; i<nSpecie_; i++)
460 const scalar Yi = Y_[i][celli];
461 c_[i] = rhoi*Yi/specieThermo_[i].W();
464 omega(c_, Ti, pi, dcdt_);
466 for (
label i=0; i<nSpecie_; i++)
468 RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
474 template<
class ReactionThermo,
class ThermoType>
475 template<
class DeltaTType>
478 const DeltaTType& deltaT
483 scalar deltaTMin = great;
485 if (!this->chemistry_)
500 scalar Ti = T[celli];
504 const scalar rhoi = rho[celli];
505 scalar
pi = p[celli];
507 for (
label i=0; i<nSpecie_; i++)
509 c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
514 scalar timeLeft = deltaT[celli];
517 while (timeLeft > small)
519 scalar dt = timeLeft;
520 this->
solve(c_, Ti, pi, dt, this->deltaTChem_[celli]);
524 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
526 this->deltaTChem_[celli] =
527 min(this->deltaTChem_[celli], this->deltaTChemMax_);
529 for (
label i=0; i<nSpecie_; i++)
532 (c_[i] - c0[i])*specieThermo_[i].
W()/deltaT[celli];
537 for (
label i=0; i<nSpecie_; i++)
548 template<
class ReactionThermo,
class ThermoType>
563 template<
class ReactionThermo,
class ThermoType>
569 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/s^3].
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
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
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...
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
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.
StandardChemistryModel(ReactionThermo &thermo)
Construct from thermo.
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(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-7/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 List< specieCoeffs > & rhs() const
Return the components of the right hand side.
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...
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