Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reorganize reaction rate classes #1223

Merged
merged 7 commits into from
Mar 20, 2022
Merged

Reorganize reaction rate classes #1223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d3522b2
[Kinetics] Make ReactionFactory internal
speth Mar 19, 2022
8d887de
[Kinetics] Put ReactionData implementations with matching rates
speth Mar 20, 2022
65aaabd
[Kinetics] Move TwoTempPlasma to its own header file
speth Mar 20, 2022
7a85750
[Kinetics] Move BlowersMaselRate to its own header file
speth Mar 20, 2022
201d556
[Kinetics] Move PlogRate to its own header file
speth Mar 20, 2022
7529467
[Kinetics] Move ChebyshevRate to its own header file
speth Mar 20, 2022
ac14f3d
[Kinetics] Deprecate obsolete reaction rate classes/methods
speth Mar 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
219 changes: 10 additions & 209 deletions include/cantera/kinetics/Arrhenius.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,16 @@ class ArrheniusBase : public ReactionRate
Units m_rate_units; //!< Reaction rate units
};

//! Data container holding shared data specific to ArrheniusRate
/**
* The data container `ArrheniusData` holds precalculated data common to
* all `ArrheniusRate` objects.
*/
struct ArrheniusData : public ReactionData
{
virtual bool update(const ThermoPhase& phase, const Kinetics& kin);
using ReactionData::update;
};

//! Arrhenius reaction rate type depends only on temperature
/*!
Expand Down Expand Up @@ -203,213 +213,6 @@ class Arrhenius3 : public ArrheniusBase
};


//! Two temperature plasma reaction rate type depends on both
//! gas temperature and electron temperature.
/*!
* The form of the two temperature plasma reaction rate coefficient is similar to an
* Arrhenius reaction rate coefficient. The temperature exponent (b) is applied to
* the electron temperature instead. In addition, the exponential term with
* activation energy for electron is included.
*
* \f[
* k_f = A T_e^b \exp (-E_a/RT) \exp (-E_{a,e}/RT_e)
* \f]
*
* Ref.: Kossyi, I. A., Kostinsky, A. Y., Matveyev, A. A., & Silakov, V. P. (1992).
* Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures.
* Plasma Sources Science and Technology, 1(3), 207.
* doi: 10.1088/0963-0252/1/3/011
*
* @ingroup arrheniusGroup
*/
class TwoTempPlasma : public ArrheniusBase
{
public:
TwoTempPlasma();

//! Constructor.
/*!
* @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units
* depend on the reaction order and the dimensionality (surface or bulk).
* @param b Temperature exponent (non-dimensional)
* @param Ea Activation energy in energy units [J/kmol]
* @param EE Activation electron energy in energy units [J/kmol]
*/
TwoTempPlasma(double A, double b, double Ea=0.0, double EE=0.0);

virtual const std::string type() const override {
return "two-temperature-plasma";
}

virtual void setContext(const Reaction& rxn, const Kinetics& kin) override;

//! Evaluate reaction rate
/*!
* @param shared_data data shared by all reactions of a given type
*/
double evalFromStruct(const TwoTempPlasmaData& shared_data) const {
// m_E4_R is the electron activation (in temperature units)
return m_A * std::exp(m_b * shared_data.logTe -
m_Ea_R * shared_data.recipT +
m_E4_R * (shared_data.electronTemp - shared_data.temperature)
* shared_data.recipTe * shared_data.recipT);
}

//! Evaluate derivative of reaction rate with respect to temperature
//! divided by reaction rate
/*!
* This method does not consider changes of electron temperature.
* A corresponding warning is raised.
* @param shared_data data shared by all reactions of a given type
*/
double ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const;

//! Return the electron activation energy *Ea* [J/kmol]
double activationElectronEnergy() const {
return m_E4_R * GasConstant;
}
};


//! Blowers Masel reaction rate type depends on the enthalpy of reaction
/**
* The Blowers Masel approximation is written by Paul Blowers,
* Rich Masel (DOI: https://doi.org/10.1002/aic.690461015) to
* adjust the activation energy based on enthalpy change of a reaction:
*
* \f{eqnarray*}{
* E_a &=& 0\; \text{if }\Delta H < -4E_0 \\
* E_a &=& \Delta H\; \text{if }\Delta H > 4E_0 \\
* E_a &=& \frac{(w + \Delta H / 2)(V_P - 2w +
* \Delta H)^2}{(V_P^2 - 4w^2 + (\Delta H)^2)}\; \text{Otherwise}
* \f}
* where
* \f[
* V_P = \frac{2w (w + E_0)}{w - E_0},
* \f]
* \f$ w \f$ is the average bond dissociation energy of the bond breaking
* and that being formed in the reaction. Since the expression is
* very insensitive to \f$ w \f$ for \f$ w >= 2 E_0 \f$, \f$ w \f$
* can be approximated to an arbitrary high value like 1000 kJ/mol.
*
* After the activation energy is determined by Blowers-Masel approximation,
* it can be plugged into Arrhenius function to calculate the rate constant.
* \f[
* k_f = A T^b \exp (-E_a/RT)
* \f]
*
* @ingroup arrheniusGroup
*/
class BlowersMasel : public ArrheniusBase
{
public:
//! Default constructor.
BlowersMasel();

//! Constructor.
/*!
* @param A Pre-exponential factor. The unit system is (kmol, m, s); actual units
* depend on the reaction order and the dimensionality (surface or bulk).
* @param b Temperature exponent (non-dimensional)
* @param Ea0 Intrinsic activation energy in energy units [J/kmol]
* @param w Average bond dissociation energy of the bond being formed and
* broken in the reaction, in energy units [J/kmol]
*/
BlowersMasel(double A, double b, double Ea0, double w);

virtual const std::string type() const override {
return "Blowers-Masel";
}

virtual void setContext(const Reaction& rxn, const Kinetics& kin) override;

//! Evaluate reaction rate
double evalRate(double logT, double recipT) const {
double Ea_R = effectiveActivationEnergy_R(m_deltaH_R);
return m_A * std::exp(m_b * logT - Ea_R * recipT);
}

//! Update information specific to reaction
void updateFromStruct(const BlowersMaselData& shared_data) {
if (shared_data.ready) {
m_deltaH_R = 0.;
for (const auto& item : m_stoich_coeffs) {
m_deltaH_R += shared_data.partialMolarEnthalpies[item.first] * item.second;
}
m_deltaH_R /= GasConstant;
}
}

//! Evaluate reaction rate
/*!
* @param shared_data data shared by all reactions of a given type
*/
double evalFromStruct(const BlowersMaselData& shared_data) const {
double Ea_R = effectiveActivationEnergy_R(m_deltaH_R);
return m_A * std::exp(m_b * shared_data.logT - Ea_R * shared_data.recipT);
}

//! divided by reaction rate
/*!
* This method does not consider potential changes due to a changed reaction
* enthalpy. A corresponding warning is raised.
* @param shared_data data shared by all reactions of a given type
*/
double ddTScaledFromStruct(const BlowersMaselData& shared_data) const;

protected:
//! Return the effective activation energy (a function of the delta H of reaction)
//! divided by the gas constant (i.e. the activation temperature) [K]
//! @internal The enthalpy change of reaction is not an independent parameter
double effectiveActivationEnergy_R(double deltaH_R) const {
if (deltaH_R < -4 * m_Ea_R) {
return 0.;
}
if (deltaH_R > 4 * m_Ea_R) {
return deltaH_R;
}
// m_E4_R is the bond dissociation energy "w" (in temperature units)
double vp = 2 * m_E4_R * ((m_E4_R + m_Ea_R) / (m_E4_R - m_Ea_R)); // in Kelvin
double vp_2w_dH = (vp - 2 * m_E4_R + deltaH_R); // (Vp - 2 w + dH)
return (m_E4_R + deltaH_R / 2) * (vp_2w_dH * vp_2w_dH) /
(vp * vp - 4 * m_E4_R * m_E4_R + deltaH_R * deltaH_R); // in Kelvin
}

public:
virtual double activationEnergy() const override {
return effectiveActivationEnergy_R(m_deltaH_R) * GasConstant;
}

//! Return the bond dissociation energy *w* [J/kmol]
double bondEnergy() const {
return m_E4_R * GasConstant;
}

//! Return current enthalpy change of reaction [J/kmol]
double deltaH() const {
return m_deltaH_R * GasConstant;
}

//! Set current enthalpy change of reaction [J/kmol]
/*!
* @internal used for testing purposes only; note that this quantity is not an
* independent variable and will be overwritten during an update of the state.
*
* @warning This method is an experimental part of the %Cantera API and
* may be changed or removed without notice.
*/
void setDeltaH(double deltaH) {
m_deltaH_R = deltaH / GasConstant;
}

protected:
//! Pairs of species indices and multiplers to calculate enthalpy change
std::vector<std::pair<size_t, double>> m_stoich_coeffs;

double m_deltaH_R; //!< enthalpy change of reaction (in temperature units)
};


//! A class template for bulk phase reaction rate specifications
template <class RateType, class DataType>
class BulkRate : public RateType
Expand Down Expand Up @@ -456,8 +259,6 @@ class BulkRate : public RateType
};

typedef BulkRate<Arrhenius3, ArrheniusData> ArrheniusRate;
typedef BulkRate<TwoTempPlasma, TwoTempPlasmaData> TwoTempPlasmaRate;
typedef BulkRate<BlowersMasel, BlowersMaselData> BlowersMaselRate;

}

Expand Down
Loading