33 template<
class ReactionThermo,
class ThermoType>
36 const ReactionThermo& thermo
45 (this->
thermo()).specieThermos()
59 nReaction_(reactions_.size()),
82 "RR." + Y_[fieldi].
name(),
94 Info<<
"StandardChemistryModel: Number of species = " << nSpecie_
95 <<
" and reactions = " << nReaction_ <<
endl;
101 template<
class ReactionThermo,
class ThermoType>
109 template<
class ReactionThermo,
class ThermoType>
126 R.
omega(p, T, c, li, dcdt);
131 template<
class ReactionThermo,
class ThermoType>
148 scalar w = R.
omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef);
153 template<
class ReactionThermo,
class ThermoType>
162 const scalar T = c[nSpecie_];
163 const scalar p = c[nSpecie_ + 1];
167 c_[i] =
max(c[i], 0);
170 omega(p, T, c_, li, dcdt);
176 for (
label i = 0; i < nSpecie_; i++)
178 const scalar
W = specieThermos_[i].W();
183 for (
label i=0; i<nSpecie_; i++)
185 cp += c_[i]*specieThermos_[i].cp(p, T);
190 for (
label i = 0; i < nSpecie_; i++)
192 const scalar hi = specieThermos_[i].ha(p, T);
197 dcdt[nSpecie_] = -dT;
200 dcdt[nSpecie_ + 1] = 0;
204 template<
class ReactionThermo,
class ThermoType>
214 const scalar T = c[nSpecie_];
215 const scalar p = c[nSpecie_ + 1];
219 c_[i] =
max(c[i], 0);
229 for (
label i = 0; i < nSpecie_; i++)
231 hi[i] = specieThermos_[i].ha(p, T);
232 cpi[i] = specieThermos_[i].cp(p, T);
240 R.
dwdc(p, T, c_, li, J, dcdt, omegaI, kfwd, kbwd,
false, dummy);
241 R.
dwdT(p, T, c_, li, omegaI, kfwd, kbwd, J,
false, dummy, nSpecie_);
247 scalar dcpdTMean = 0;
248 for (
label i=0; i<nSpecie_; i++)
250 cpMean += c_[i]*cpi[i];
251 dcpdTMean += c_[i]*specieThermos_[i].dcpdT(p, T);
254 for (
label i=0; i<nSpecie_; i++)
256 dTdt += hi[i]*dcdt[i];
260 for (
label i = 0; i < nSpecie_; i++)
263 for (
label j = 0; j < nSpecie_; j++)
265 J(nSpecie_, i) += hi[j]*J(j, i);
267 J(nSpecie_, i) += cpi[i]*dTdt;
268 J(nSpecie_, i) /= -cpMean;
272 J(nSpecie_, nSpecie_) = 0;
273 for (
label i = 0; i < nSpecie_; i++)
275 J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
277 J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
278 J(nSpecie_, nSpecie_) /= -cpMean;
279 J(nSpecie_, nSpecie_) += dTdt/
T;
283 template<
class ReactionThermo,
class ThermoType>
294 extrapolatedCalculatedFvPatchScalarField::typeName
306 const label nReaction = reactions_.size();
308 scalar pf, cf, pr, cr;
311 if (this->chemistry_)
315 const scalar rhoi = rho[celli];
316 const scalar Ti = T[celli];
317 const scalar
pi = p[celli];
321 for (
label i=0; i<nSpecie_; i++)
323 c_[i] = rhoi*Y_[i][celli]/specieThermos_[i].W();
331 R.
omega(pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef);
335 tc[celli] += R.
rhs()[
s].stoichCoeff*pf*cf;
339 tc[celli] = nReaction*cSum/tc[celli];
343 ttc.
ref().correctBoundaryConditions();
349 template<
class ReactionThermo,
class ThermoType>
363 if (this->chemistry_)
371 const scalar hi = specieThermos_[i].Hf();
372 Qdot[celli] -= hi*RR_[i][celli];
381 template<
class ReactionThermo,
class ThermoType>
407 scalar pf, cf, pr, cr;
412 const scalar rhoi = rho[celli];
413 const scalar Ti = T[celli];
414 const scalar
pi = p[celli];
416 for (
label i=0; i<nSpecie_; i++)
418 const scalar Yi = Y_[i][celli];
419 c_[i] = rhoi*Yi/specieThermos_[i].W();
423 const scalar omegai = R.
omega 425 pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef
430 if (si == R.
lhs()[
s].index)
432 RR[celli] -= R.
lhs()[
s].stoichCoeff*omegai;
438 if (si == R.
rhs()[
s].index)
440 RR[celli] += R.
rhs()[
s].stoichCoeff*omegai;
444 RR[celli] *= specieThermos_[si].W();
451 template<
class ReactionThermo,
class ThermoType>
454 if (!this->chemistry_)
467 const scalar rhoi = rho[celli];
468 const scalar Ti = T[celli];
469 const scalar
pi = p[celli];
471 for (
label i=0; i<nSpecie_; i++)
473 const scalar Yi = Y_[i][celli];
474 c_[i] = rhoi*Yi/specieThermos_[i].W();
477 omega(pi, Ti, c_, celli, dcdt_);
479 for (
label i=0; i<nSpecie_; i++)
481 RR_[i][celli] = dcdt_[i]*specieThermos_[i].W();
487 template<
class ReactionThermo,
class ThermoType>
488 template<
class DeltaTType>
491 const DeltaTType& deltaT
496 scalar deltaTMin = great;
498 if (!this->chemistry_)
513 scalar Ti = T[celli];
517 const scalar rhoi = rho[celli];
518 scalar
pi = p[celli];
520 for (
label i=0; i<nSpecie_; i++)
522 c_[i] = rhoi*Y_[i][celli]/specieThermos_[i].W();
527 scalar timeLeft = deltaT[celli];
530 while (timeLeft > small)
532 scalar dt = timeLeft;
533 this->
solve(pi, Ti, c_, celli, dt, this->deltaTChem_[celli]);
537 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
539 this->deltaTChem_[celli] =
540 min(this->deltaTChem_[celli], this->deltaTChemMax_);
542 for (
label i=0; i<nSpecie_; i++)
545 (c_[i] - c0[i])*specieThermos_[i].
W()/deltaT[celli];
550 for (
label i=0; i<nSpecie_; i++)
561 template<
class ReactionThermo,
class ThermoType>
576 template<
class ReactionThermo,
class ThermoType>
582 return this->solve<scalarField>(deltaT);
StandardChemistryModel(const ReactionThermo &thermo)
Construct from thermo.
#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 > &)
basicSpecieMixture & composition
T & ref() const
Return non-const reference or generate a fatal error.
void dwdc(const scalar p, const scalar T, const scalarField &c, const label li, 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 void jacobian(const scalar t, const scalarField &c, const label li, scalarField &dcdt, scalarSquareMatrix &J) const
Calculate the Jacobian of the system.
virtual void omega(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dcdt) const
dc/dt = omega, rate of change in concentration, for each species
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool cp(const fileName &src, const fileName &dst, const bool followLink=true)
Copy, recursively if necessary, the source to the destination.
tmp< volScalarField > trho
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
rhoReactionThermo & thermo
Extends base chemistry model by adding a thermo package, and ODE functions. Introduces chemistry equa...
void omega(const scalar p, const scalar T, const scalarField &c, const label li, scalarField &dcdt) const
Net reaction rate for individual species.
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.
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
Foam::multiComponentMixture.
virtual scalar omegaI(label iReaction, const scalar p, const scalar T, const scalarField &c, const label li, scalar &pf, scalar &cf, label &lRef, scalar &pr, scalar &cr, label &rRef) const
Return the reaction rate for iReaction and the reference.
virtual ~StandardChemistryModel()
Destructor.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-8/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 void derivatives(const scalar t, const scalarField &c, const label li, scalarField &dcdt) const
Calculate the derivatives in dydx.
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)
void dwdT(const scalar p, const scalar T, const scalarField &c, const label li, 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.
A class for managing temporary objects.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...