32 template<
class CompType,
class Sol
idThermo,
class GasThermo>
36 typename CompType::reactionThermo& thermo
47 gasThermo_(pyrolisisGases_.size()),
48 nGases_(pyrolisisGases_.size()),
49 nSpecie_(this->Ys_.size() + nGases_),
59 this->Ys_[fieldi].
name() +
"0",
75 this->Ys_[fieldi].
name() +
"0",
107 this->Ys_[fieldi].
name() +
"0",
118 Ys0_[fieldi].primitiveFieldRef() =
120 *
max(this->Ys_[fieldi], scalar(0.001))*this->
mesh().V();
133 "RRg." + pyrolisisGases_[fieldi],
148 this->
mesh().template lookupObject<dictionary>
151 ).subDict(pyrolisisGases_[gasI]);
156 new GasThermo(thermoDict)
160 Info<<
"pyrolysisChemistryModel: " <<
nl;
161 Info<<
indent <<
"Number of solids = " << this->nSolids_ <<
nl;
163 forAll(this->reactions_, i)
165 Info<< dynamic_cast<const SolidReaction<SolidThermo>&>
175 template<
class CompType,
class Sol
idThermo,
class GasThermo>
183 template<
class CompType,
class Sol
idThermo,
class GasThermo>
193 scalar pf, cf, pr, cr;
196 const label celli = cellCounter_;
200 forAll(this->reactions_, i)
208 scalar omegai = omega
210 R, c, T, p, pf, cf, lRef, pr, cr, rRef
217 rhoL = this->solidThermo_[si].rho(p, T);
223 scalar rhoR = this->solidThermo_[si].rho(p, T);
229 Ys0_[si][celli] += sr*omegai;
235 om[gi + this->nSolids_] += (1.0 - sr)*omegai;
243 template<
class CompType,
class Sol
idThermo,
class GasThermo>
261 label celli = cellCounter_;
263 for (
label i=0; i<nSpecie_; i++)
265 c1[i] =
max(0.0, c[i]);
268 scalar kf = R.
kf(p, T, c1);
275 const scalar
exp = R.
lhs()[si].exponent;
278 pow(c1[si]/Ys0_[si][celli], exp)
286 template<
class CompType,
class Sol
idThermo,
class GasThermo>
304 scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
309 template<
class CompType,
class Sol
idThermo,
class GasThermo>
318 const scalar T = c[nSpecie_];
319 const scalar p = c[nSpecie_ + 1];
323 dcdt = omega(c, T, p);
327 for (
label i=0; i<this->nSolids_; i++)
334 for (
label i=0; i<this->nSolids_; i++)
336 scalar dYidt = dcdt[i]/cTot;
337 scalar Yi = c[i]/cTot;
338 newCp += Yi*this->solidThermo_[i].Cp(p, T);
339 newhi -= dYidt*this->solidThermo_[i].Hc();
342 scalar dTdt = newhi/newCp;
343 scalar dtMag =
min(500.0,
mag(dTdt));
344 dcdt[nSpecie_] = dTdt*dtMag/(
mag(dTdt) + 1.0e-10);
347 dcdt[nSpecie_ + 1] = 0.0;
351 template<
class CompType,
class Sol
idThermo,
class GasThermo>
361 const scalar T = c[nSpecie_];
362 const scalar p = c[nSpecie_ + 1];
366 for (
label i=0; i<this->nSolids_; i++)
368 c2[i] =
max(c[i], 0.0);
371 for (
label i=0; i<nEqns(); i++)
373 for (
label j=0; j<nEqns(); j++)
380 dcdt = omega(c2, T, p);
382 for (
label ri=0; ri<this->reactions_.size(); ri++)
386 scalar kf0 = R.
kf(p, T, c2);
395 scalar
exp = R.
lhs()[i].exponent;
402 kf *= exp*
pow(c2[si], exp - 1.0);
411 kf *= exp*
pow(c2[si], exp - 1.0);
416 Info<<
"Solid reactions have only elements on slhs" 436 scalar
delta = 1.0e-8;
440 for (
label i=0; i<nEqns(); i++)
442 dfdc[i][nSpecie_] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
448 template<
class CompType,
class Sol
idThermo,
class GasThermo>
453 return (nSpecie_ + 2);
457 template<
class CompType,
class Sol
idThermo,
class GasThermo>
461 if (!this->chemistry_)
482 this->RRs_[i].field() = 0.0;
487 RRg_[i].field() = 0.0;
492 cellCounter_ = celli;
494 const scalar
delta = this->
mesh().V()[celli];
496 if (this->reactingCells_[celli])
498 scalar rhoi = rho[celli];
503 for (
label i=0; i<this->nSolids_; i++)
505 c[i] = rhoi*this->Ys_[i][celli]*
delta;
512 this->RRs_[i][celli] = dcdt[i]/
delta;
517 RRg_[i][celli] = dcdt[this->nSolids_ + i]/
delta;
524 template<
class CompType,
class Sol
idThermo,
class GasThermo>
531 scalar deltaTMin = great;
533 if (!this->chemistry_)
554 this->RRs_[i].field() = 0.0;
558 RRg_[i].field() = 0.0;
571 if (this->reactingCells_[celli])
573 cellCounter_ = celli;
575 scalar rhoi = rho[celli];
576 scalar
pi = p[celli];
577 scalar Ti = T[celli];
579 for (
label i=0; i<this->nSolids_; i++)
581 c[i] = rhoi*this->Ys_[i][celli]*delta[celli];
587 scalar timeLeft = deltaT;
590 while (timeLeft > small)
592 scalar dt = timeLeft;
593 this->
solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
597 deltaTMin =
min(this->deltaTChem_[celli], deltaTMin);
602 this->RRs_[i][celli] = dc[i]/(deltaT*delta[celli]);
607 RRg_[i][celli] = dc[this->nSolids_ + i]/(deltaT*delta[celli]);
611 dc = omega(c0, Ti, pi,
true);
616 deltaTMin =
min(deltaTMin, 2*deltaT);
622 template<
class CompType,
class Sol
idThermo,
class GasThermo>
637 "Hs_" + pyrolisisGases_[index],
638 this->mesh_.time().timeName(),
651 const GasThermo&
mixture = gasThermo_[index];
655 gasHs[celli] = mixture.Hs(p[celli], T[celli]);
662 template<
class CompType,
class Sol
idThermo,
class GasThermo>
virtual ~pyrolysisChemistryModel()
Destructor.
virtual const speciesTable & gasSpecies() const
Access to gas specie list.
#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.
bool typeHeaderOk(const bool checkType=true)
Read header (uses typeFilePath to find file) and check header.
Ostream & indent(Ostream &os)
Indent stream.
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
A list of keyword definitions, which are a keyword followed by any number of values (e...
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.
T & ref() const
Return non-const reference or generate a fatal error.
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
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.
virtual const volScalarField & T() const
Temperature [K].
rhoReactionThermo & thermo
Extends base solid chemistry model by adding a thermo package, and ODE functions. Introduces chemistr...
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
virtual volScalarField & p()
Pressure [Pa].
virtual const List< specieCoeffs > & grhs() const
Access to the gas components of the left hand side.
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.
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))
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].
virtual scalar solve(const scalar deltaT)
Solve the reaction system for the given time step.
const word dictName("particleTrackDict")
Fundamental solid thermodynamic properties.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
const dimensionSet dimEnergy
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
pyrolysisChemistryModel(typename CompType::reactionThermo &thermo)
Construct from thermo.
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...
Read solid reactions of the type S1 = S2 + G1.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
const dimensionedVector & g
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
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.
virtual scalar kf(const scalar p, const scalar T, const scalarField &c) const =0
Forward rate constant.
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.