35 template<
class CompType,
class Sol
idThermo,
class GasThermo>
39 typename CompType::reactionThermo&
thermo 43 pyrolisisGases_(this->reactions_[0].gasSpecies()),
44 gasThermo_(pyrolisisGases_.size()),
45 nGases_(pyrolisisGases_.size()),
46 nSpecie_(this->Ys_.size() + nGases_),
104 this->
Ys_[fieldi].
name() +
"0",
115 Ys0_[fieldi].primitiveFieldRef() =
117 *
max(this->
Ys_[fieldi], scalar(0.001))*this->
mesh().V();
145 this->
mesh().template lookupObject<dictionary>
157 Info<<
"pyrolysisChemistryModel: " <<
nl;
162 Info<< dynamic_cast<const solidReaction<SolidThermo>&>
172 template<
class CompType,
class Sol
idThermo,
class GasThermo>
180 template<
class CompType,
class Sol
idThermo,
class GasThermo>
190 scalar pf, cf, pr, cr;
193 const label celli = cellCounter_;
197 forAll(this->reactions_, i)
201 scalar omegai = omega
203 R,
c,
T,
p, pf, cf, lRef, pr, cr, rRef
208 label si =
R.lhs()[
s].index;
210 rhoL = this->solidThermo_[si].rho(
p,
T);
215 label si =
R.rhs()[
s].index;
216 scalar rhoR = this->solidThermo_[si].rho(
p,
T);
222 Ys0_[si][celli] += sr*omegai;
227 label gi =
R.grhs()[
g].index;
228 om[gi + this->nSolids_] += (1.0 - sr)*omegai;
236 template<
class CompType,
class Sol
idThermo,
class GasThermo>
254 label celli = cellCounter_;
256 for (label i=0; i<nSpecie_; i++)
261 scalar kf =
R.kf(
p,
T,
c1);
263 const label Nl =
R.lhs().size();
265 for (label
s=0;
s<Nl;
s++)
267 label si =
R.lhs()[
s].index;
268 const scalar
exp =
R.lhs()[
s].exponent;
279 template<
class CompType,
class Sol
idThermo,
class GasThermo>
297 scalar w = omega(
R,
c,
T,
p, pf, cf, lRef, pr, cr, rRef);
302 template<
class CompType,
class Sol
idThermo,
class GasThermo>
311 const scalar
T =
c[nSpecie_];
312 const scalar
p =
c[nSpecie_ + 1];
316 dcdt = omega(
c,
T,
p);
320 for (label i=0; i<this->nSolids_; i++)
327 for (label i=0; i<this->nSolids_; i++)
329 scalar dYidt = dcdt[i]/cTot;
330 scalar Yi =
c[i]/cTot;
331 newCp += Yi*this->solidThermo_[i].Cp(
p,
T);
332 newhi -= dYidt*this->solidThermo_[i].Hc();
335 scalar dTdt = newhi/newCp;
336 scalar dtMag =
min(500.0,
mag(dTdt));
337 dcdt[nSpecie_] = dTdt*dtMag/(
mag(dTdt) + 1.0e-10);
340 dcdt[nSpecie_ + 1] = 0.0;
344 template<
class CompType,
class Sol
idThermo,
class GasThermo>
354 const scalar
T =
c[nSpecie_];
355 const scalar
p =
c[nSpecie_ + 1];
359 for (label i=0; i<this->nSolids_; i++)
364 for (label i=0; i<nEqns(); i++)
366 for (label j=0; j<nEqns(); j++)
373 dcdt = omega(
c2,
T,
p);
375 for (label ri=0; ri<this->reactions_.size(); ri++)
377 const Reaction<SolidThermo>&
R = this->reactions_[ri];
379 scalar kf0 =
R.kf(
p,
T,
c2);
383 label sj =
R.lhs()[j].index;
387 label si =
R.lhs()[i].index;
388 scalar
exp =
R.lhs()[i].exponent;
409 Info<<
"Solid reactions have only elements on slhs" 417 label si =
R.lhs()[i].index;
422 label si =
R.rhs()[i].index;
429 scalar
delta = 1.0e-8;
433 for (label i=0; i<nEqns(); i++)
435 dfdc[i][nSpecie_] = 0.5*(dcdT1[i] - dcdT0[i])/
delta;
441 template<
class CompType,
class Sol
idThermo,
class GasThermo>
446 return (nSpecie_ + 2);
450 template<
class CompType,
class Sol
idThermo,
class GasThermo>
454 if (!this->chemistry_)
470 this->solidThermo().
rho()
475 this->RRs_[i].field() = 0.0;
480 RRg_[i].field() = 0.0;
485 cellCounter_ = celli;
489 if (this->reactingCells_[celli])
491 scalar rhoi =
rho[celli];
492 scalar Ti = this->solidThermo().T()[celli];
493 scalar
pi = this->solidThermo().p()[celli];
496 for (label i=0; i<this->nSolids_; i++)
498 c[i] = rhoi*this->Ys_[i][celli]*
delta;
505 this->RRs_[i][celli] = dcdt[i]/
delta;
510 RRg_[i][celli] = dcdt[this->nSolids_ + i]/
delta;
517 template<
class CompType,
class Sol
idThermo,
class GasThermo>
524 scalar deltaTMin = GREAT;
526 if (!this->chemistry_)
542 this->solidThermo().
rho()
547 this->RRs_[i].field() = 0.0;
551 RRg_[i].field() = 0.0;
564 if (this->reactingCells_[celli])
566 cellCounter_ = celli;
568 scalar rhoi =
rho[celli];
569 scalar
pi =
p[celli];
570 scalar Ti =
T[celli];
572 for (label i=0; i<this->nSolids_; i++)
574 c[i] = rhoi*this->Ys_[i][celli]*
delta[celli];
580 scalar timeLeft = deltaT;
583 while (timeLeft > SMALL)
585 scalar dt = timeLeft;
586 this->
solve(c, Ti,
pi, dt, this->deltaTChem_[celli]);
590 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
595 this->RRs_[i][celli] = dc[i]/(deltaT*
delta[celli]);
600 RRg_[i][celli] = dc[this->nSolids_ + i]/(deltaT*
delta[celli]);
604 dc = omega(c0, Ti,
pi,
true);
609 deltaTMin =
min(deltaTMin, 2*deltaT);
615 template<
class CompType,
class Sol
idThermo,
class GasThermo>
626 "Hs_" + pyrolisisGases_[index],
631 auto& gasHs = tHs.ref();
633 const GasThermo&
mixture = gasThermo_[index];
637 gasHs[celli] =
mixture.Hs(
p[celli],
T[celli]);
644 template<
class CompType,
class Sol
idThermo,
class GasThermo>
word dictName() const
The local dictionary name (final part of scoped name)
virtual ~pyrolysisChemistryModel()
Destructor.
Ostream & indent(Ostream &os)
Indent stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
PtrList< GasThermo > gasThermo_
Thermodynamic data of gases.
A list of keyword definitions, which are a keyword followed by a number of values (eg...
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.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
constexpr char nl
The newline '\n' character (0x0a)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
speciesTable pyrolisisGases_
List of gas species present in reaction system.
Ignore writing from objectRegistry::writeObject()
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
Extends base solid chemistry model by adding a thermo package, and ODE functions. ...
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
PtrList< volScalarField::Internal > RRs_
List of reaction rate per solid [kg/m3/s].
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
virtual scalarField omega(const scalarField &c, const scalar T, const scalar p, const bool updateC0=false) const
dc/dt = omega, rate of change in concentration, for each species
const dimensionSet dimVolume(pow3(dimLength))
virtual label nEqns() const
Number of ODE's to solve.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
dimensionedScalar exp(const dimensionedScalar &ds)
virtual tmp< volScalarField > gasHs(const volScalarField &p, const volScalarField &T, const label i) const
Return sensible enthalpy for gas i [J/Kg].
label nGases_
Number of gas species.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual scalar solve(const scalar deltaT)
Solve the reaction system for the given time step.
label nSolids_
Number of solid components.
constexpr scalar pi(M_PI)
const dictionary & thermoDict
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
PtrList< volScalarField::Internal > RRg_
List of reaction rate per gas [kg/m3/s].
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
const PtrList< Reaction< SolidThermo > > & reactions_
Reactions.
SolverPerformance< Type > solve(faMatrix< Type > &, const dictionary &solverControls)
Solve returning the solution statistics given convergence tolerance.
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (respects is_globalIOobject trait) and check its info.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
const uniformDimensionedVectorField & g
PtrList< volScalarField > Ys0_
List of accumulative solid concentrations.
Fundamental solid thermodynamic properties.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const dimensionSet dimEnergy
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
pyrolysisChemistryModel(typename CompType::reactionThermo &thermo)
Construct from thermo.
PtrList< volScalarField > & Ys_
Reference to solid mass fractions.
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
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...
Automatically write from objectRegistry::writeObject()
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
Information stream (stdout output on master, null elsewhere)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Defines the attributes of an object for which implicit objectRegistry management is supported...
virtual void calculate()
Calculates the reaction rates.
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
Do not request registration (bool: false)
static constexpr const zero Zero
Global zero (0)