13 #include <boost/math/tools/roots.hpp> 16 namespace bmt = boost::math::tools;
30 void Reactor::setKineticsMgr(Kinetics& kin)
44 "Error: reactor is empty.");
69 for (
auto& S : m_surfaces) {
70 S->getCoverages(y + loc);
71 loc += S->thermo()->nSpecies();
77 if (!m_thermo || (m_chem && !
m_kin)) {
78 throw CanteraError(
"Reactor::initialize",
"Reactor contents not set" 79 " for reactor '" + m_name +
"'.");
89 for (
size_t n = 0; n < m_wall.size(); n++) {
96 for (
auto& S : m_surfaces) {
97 m_nv += S->thermo()->nSpecies();
98 size_t nt = S->kinetics()->nTotalSpecies();
99 maxnt = std::max(maxnt, nt);
100 if (m_chem && &
m_kin->
thermo(0) != &S->kinetics()->thermo(0)) {
102 "First phase of all kinetics managers must be the gas.");
105 m_work.resize(maxnt);
110 size_t ns = m_sensParams.size();
111 for (
auto& S : m_surfaces) {
112 ns += S->nSensParams();
135 auto u_err = [
this, U](
double T) {
141 boost::uintmax_t maxiter = 100;
142 std::pair<double, double> TT;
144 TT = bmt::bracket_and_solve_root(
145 u_err, T, 1.2,
true, bmt::eps_tolerance<double>(48), maxiter);
146 }
catch (std::exception& err) {
151 bmt::eps_tolerance<double>(48), maxiter);
152 }
catch (std::exception& err2) {
156 "{}\nat U = {}, rho = {}", err2.what(), U,
m_mass / m_vol);
159 if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
160 throw CanteraError(
"Reactor::updateState",
"root finding failed");
179 for (
auto& S : m_surfaces) {
180 S->setCoverages(y+loc);
181 loc += S->thermo()->nSpecies();
186 doublereal* ydot, doublereal* params)
189 double* dYdt = ydot + 3;
207 for (
size_t k = 0; k <
m_nsp; k++) {
211 dYdt[k] -= Y[k] * mdot_surf /
m_mass;
225 for (
size_t i = 0; i < m_outlet.size(); i++) {
226 double mdot_out = m_outlet[i]->massFlowRate(time);
229 ydot[2] -= mdot_out * m_enthalpy;
234 for (
size_t i = 0; i < m_inlet.size(); i++) {
235 double mdot_in = m_inlet[i]->massFlowRate(time);
237 for (
size_t n = 0; n <
m_nsp; n++) {
238 double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
240 dYdt[n] += (mdot_spec - mdot_in * Y[n]) /
m_mass;
243 ydot[2] += mdot_in * m_inlet[i]->enthalpy_mass();
255 for (
size_t i = 0; i < m_wall.size(); i++) {
256 int lr = 1 - 2*m_lr[i];
257 m_vdot += lr*m_wall[i]->vdot(t);
258 m_Q += lr*m_wall[i]->Q(t);
267 double mdot_surf = 0.0;
269 for (
auto S : m_surfaces) {
281 for (
size_t k = 1; k < nk; k++) {
282 ydot[loc + k] = m_work[surfloc+k]*rs0*surf->
size(k);
283 sum -= ydot[loc + k];
288 double wallarea = S->area();
289 for (
size_t k = 0; k <
m_nsp; k++) {
290 m_sdot[k] += m_work[k]*wallarea;
291 mdot_surf +=
m_sdot[k] * mw[k];
301 "Reaction number out of range ({})", rxn);
306 m_sensParams.emplace_back(
307 SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
313 throw CanteraError(
"Reactor::addSensitivitySpeciesEnthalpy",
314 "Species index out of range ({})", k);
320 m_sensParams.emplace_back(
321 SensitivityParameter{k, p, m_thermo->
Hf298SS(k),
322 SensParameterType::enthalpy});
334 size_t offset =
m_nsp;
335 for (
auto& S : m_surfaces) {
352 }
else if (nm ==
"mass") {
354 }
else if (nm ==
"volume") {
356 }
else if (nm ==
"int_energy") {
370 }
else if (k >= 3 && k <
neq()) {
372 if (k < m_thermo->nSpecies()) {
377 for (
auto& S : m_surfaces) {
379 if (k < th->nSpecies()) {
386 throw CanteraError(
"Reactor::componentName",
"Index is out of bounds.");
394 for (
auto& p : m_sensParams) {
395 if (p.type == SensParameterType::reaction) {
398 }
else if (p.type == SensParameterType::enthalpy) {
402 for (
auto& S : m_surfaces) {
403 S->setSensitivityParameters(params);
407 m_kin->invalidateCache();
416 for (
auto& p : m_sensParams) {
417 if (p.type == SensParameterType::reaction) {
419 }
else if (p.type == SensParameterType::enthalpy) {
423 for (
auto& S : m_surfaces) {
424 S->resetSensitivityParameters();
428 m_kin->invalidateCache();
virtual double size(size_t k) const
Returns the number of sites occupied by one molecule of species k.
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
void setChemistry(bool cflag=true)
Enable or disable changes in reactor composition due to chemical reactions.
vector_fp m_sdot
Production rates of gas phase species on surfaces [kmol/s].
void restoreState(const vector_fp &state)
Restore a state saved on a previous call to saveState.
doublereal m_mass
total mass
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase, assuming an ideal solution model (see Thermodynamic Properties and class SurfPhase).
virtual size_t nSensParams()
Number of sensitivity parameters associated with this reactor (including walls)
void getMassFractions(doublereal *const y) const
Get the species mass fractions.
Header file for class Wall.
vector_fp m_wdot
Species net molar production rates.
doublereal temperature() const
Temperature (K).
const doublereal * massFractions() const
Return a const pointer to the mass fraction array.
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
void saveState(vector_fp &state) const
Save the current internal state of the phase.
const size_t npos
index returned by functions to indicate "no position"
virtual void addSensitivityReaction(size_t rxn)
Add a sensitivity parameter associated with the reaction number rxn (in the homogeneous phase)...
ReactorNet & network()
The ReactorNet that this reactor belongs to.
virtual void evalWalls(double t)
Evaluate terms related to Walls.
size_t nSpecies() const
Returns the number of species in the phase.
virtual void resetSensitivity(double *params)
Reset the reaction rate multipliers.
virtual doublereal density() const
Density (kg/m^3).
virtual doublereal minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid...
doublereal enthalpy_mass() const
Specific enthalpy. Units: J/kg.
virtual size_t speciesIndex(const std::string &nm) const
Return the index in the solution vector for this reactor of the species named nm, in either the homog...
virtual std::string componentName(size_t k)
Return the name of the solution component with index i.
virtual void initialize(doublereal t0=0.0)
Initialize the reactor.
size_t registerSensitivityParameter(const std::string &name, double value, double scale)
Used by Reactor and Wall objects to register the addition of sensitivity parameters so that the React...
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Base class for a phase with thermodynamic properties.
virtual void updateState(doublereal *y)
Set the state of the reactor to correspond to the state vector y.
doublereal Hf298SS(const size_t k) const
Report the 298 K Heat of Formation of the standard state of one species (J kmol-1) ...
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected...
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
virtual void applySensitivity(double *params)
Set reaction rate multipliers based on the sensitivity variables in params.
doublereal multiplier(size_t i) const
The current value of the multiplier for reaction i.
std::string speciesName(size_t k) const
Name of the species with index k.
Represents a wall between between two ReactorBase objects.
doublereal m_Q
net heat transfer through walls [W]
virtual void syncState()
Set the state of the reactor to correspond to the state of the associated ThermoPhase object...
virtual double evalSurfaces(double t, double *ydot)
Evaluate terms related to surface reactions.
Public interface for kinetics managers.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
virtual void getState(doublereal *y)
Get the the current state of the reactor.
Base class for exceptions thrown by Cantera classes.
size_t nReactions() const
Number of reactions in the reaction mechanism.
std::string reactionString(size_t i) const
Return a string representing the reaction.
virtual void initialize()
Called just before the start of integration.
virtual void resetHf298(const size_t k=npos)
Restore the original heat of formation of one or more species.
Header file for class ReactorSurface.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
virtual void syncState()
Set the state of the reactor to correspond to the state of the associated ThermoPhase object...
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
virtual size_t componentIndex(const std::string &nm) const
Return the index in the solution vector for this reactor of the component named nm.
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
virtual void modifyOneHf298SS(const size_t k, const doublereal Hf298New)
Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1) ...
size_t m_nsp
Number of homogeneous species in the mixture.
virtual size_t neq()
Number of equations (state variables) for this reactor.
virtual void setMultiplier(size_t i, doublereal f)
Set the multiplier for reaction i to f.
virtual void evalEqs(doublereal t, doublereal *y, doublereal *ydot, doublereal *params)
doublereal m_vdot
net rate of volume change from moving walls [m^3/s]
Namespace for the Cantera kernel.
virtual void addSensitivitySpeciesEnthalpy(size_t k)
Add a sensitivity parameter associated with the enthalpy formation of species k (in the homogeneous p...
doublereal intEnergy_mass() const
Specific internal energy. Units: J/kg.
virtual doublereal maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
void setState_TR(doublereal t, doublereal rho)
Set the internally stored temperature (K) and density (kg/m^3)
virtual void updateSurfaceState(double *y)
Update the state of SurfPhase objects attached to this reactor.
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
doublereal siteDensity()
Returns the site density.
std::string name() const
Return the name of this reactor.