35template<
class CompType,
class Sol
idThermo,
class GasThermo>
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>&>
172template<
class CompType,
class Sol
idThermo,
class GasThermo>
180template<
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;
236template<
class CompType,
class Sol
idThermo,
class GasThermo>
254 label celli = cellCounter_;
256 for (label i=0; i<nSpecie_; i++)
258 c1[i] =
max(0.0, c[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;
271 pow(c1[si]/Ys0_[si][celli],
exp)
279template<
class CompType,
class Sol
idThermo,
class GasThermo>
297 scalar w = omega(
R, c,
T,
p, pf, cf, lRef, pr, cr, rRef);
302template<
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;
344template<
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++)
361 c2[i] =
max(c[i], 0.0);
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++)
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;
441template<
class CompType,
class Sol
idThermo,
class GasThermo>
446 return (nSpecie_ + 2);
450template<
class CompType,
class Sol
idThermo,
class GasThermo>
454 if (!this->chemistry_)
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];
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;
517template<
class CompType,
class Sol
idThermo,
class GasThermo>
524 scalar deltaTMin = GREAT;
526 if (!this->chemistry_)
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);
615template<
class CompType,
class Sol
idThermo,
class GasThermo>
630 "Hs_" + pyrolisisGases_[index],
631 this->mesh_.time().timeName(),
644 const GasThermo&
mixture = gasThermo_[index];
648 gasHs[celli] =
mixture.Hs(
p[celli],
T[celli]);
655template<
class CompType,
class Sol
idThermo,
class GasThermo>
#define R(A, B, C, D, E, F, K, M)
const uniformDimensionedVectorField & g
ReactionThermo reactionThermo
Thermo type.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
const T * set(const label i) const
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
virtual const volScalarField & T() const
Temperature [K].
virtual volScalarField & p()
Pressure [Pa].
static const word dictName
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
Pyrolysis chemistry model. It includes gas phase in the solid reaction.
PtrList< GasThermo > gasThermo_
Thermodynamic data of gases.
virtual tmp< volScalarField > gasHs(const volScalarField &p, const volScalarField &T, const label i) const
Return sensible enthalpy for gas i [J/Kg].
virtual void derivatives(const scalar t, const scalarField &c, scalarField &dcdt) const
Calculate the derivatives in dydx.
label nGases_
Number of gas species.
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.
virtual label nEqns() const
Number of ODE's to solve.
speciesTable pyrolisisGases_
List of gas species present in reaction system.
virtual ~pyrolysisChemistryModel()
Destructor.
PtrList< volScalarField > Ys0_
List of accumulative solid concentrations.
PtrList< volScalarField::Internal > RRg_
List of reaction rate per gas [kg/m3/s].
virtual void jacobian(const scalar t, const scalarField &c, scalarField &dcdt, scalarSquareMatrix &dfdc) const
Calculate the Jacobian of the system.
virtual void calculate()
Calculates the reaction rates.
const volScalarField & omega() const
Access functions.
Extends base solid chemistry model by adding a thermo package, and ODE functions.
PtrList< volScalarField > & Ys_
Reference to solid mass fractions.
const PtrList< Reaction< SolidThermo > > & reactions_
Reactions.
label nSolids_
Number of solid components.
PtrList< volScalarField::Internal > RRs_
List of reaction rate per solid [kg/m3/s].
Fundamental solid thermodynamic properties.
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
A class for managing temporary objects.
const dictionary & thermoDict
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
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))
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimEnergy
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ostream & indent(Ostream &os)
Indent stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVolume(pow3(dimLength))
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.