diff --git a/doc/doxygen/cantera.bib b/doc/doxygen/cantera.bib index 6e9760f8fa..d983b6b1de 100644 --- a/doc/doxygen/cantera.bib +++ b/doc/doxygen/cantera.bib @@ -31,6 +31,22 @@ @article{blowers2000 url = {https://dx.doi.org/10.1002/aic.690461015}, volume = {46}, year = {2000}} +@book{capitelli2013, + title={Plasma kinetics in atmospheric gases}, + author={Capitelli, Mario and Ferreira, Carlos M and Gordiets, Boris F and Osipov, Alexey I}, + volume={31}, + year={2013}, + publisher={Springer Science \& Business Media} +} +@article{castela2016, + title={Modelling the impact of non-equilibrium discharges on reactive mixtures for simulations of plasma-assisted ignition in turbulent flows}, + author={Castela, Maria and Fiorina, Beno{\^\i}t and Coussement, Axel and Gicquel, Olivier and Darabiha, Nasser and Laux, Christophe O}, + journal={Combustion and flame}, + volume={166}, + pages={133--147}, + year={2016}, + publisher={Elsevier} +} @book{chapman1970, title={The mathematical theory of non-uniform gases}, author={S.~Chapman and T.~G.~Cowling}, @@ -111,6 +127,16 @@ @misc{gou2011 howpublished = {\url{https://engine.princeton.edu/model-reduction/}}, note = {Accessed: 2022-04-25}, year = {2011}} +@article{guerra2019, + title={Modelling N2--O2 plasmas: volume and surface kinetics}, + author={Guerra, Vasco and Tejero-del-Caz, Antonio and Pintassilgo, Carlos D and Alves, Lu{\'\i}s L}, + journal={Plasma Sources Science and Technology}, + volume={28}, + number={7}, + pages={073001}, + year={2019}, + publisher={IOP Publishing} +} @article{gudmundsson2001, author = {J.~T.~Gudmundsson}, journal = {Plasma Sources Science and Technology}, @@ -390,6 +416,16 @@ @book{poling2001 isbn = {0-07-149999-7}, year = {2001}, edition = {Fifth}} +@article{popov2016, + title={Kinetics of plasma-assisted combustion: effect of non-equilibrium excitation on the ignition and oxidation of combustible mixtures}, + author={Popov, NA}, + journal={Plasma Sources Science and Technology}, + volume={25}, + number={4}, + pages={043002}, + year={2016}, + publisher={IOP Publishing} +} @phdthesis{prager2005, author = {J.~Prager}, school = {Technische Universität Darmstadt}, @@ -470,6 +506,16 @@ @article{singal2024 number = {1--4}, pages = {105779}, year = {2024}} +@article{starikovskiy2013, +title = "Plasma-assisted ignition and combustion", +journal = "Progress in Energy and Combustion Science", +volume = "39", +number = "1", +pages = "61 - 110", +year = "2013", +issn = "0360-1285", +author = "Andrey Starikovskiy and Nickolay Aleksandrov", +} @article{stewart1989, author = {P.~H.~Stewart and C.~W.~Larson and D.~Golden}, journal = {Combustion and Flame}, @@ -570,6 +616,19 @@ @article{westbrook1981 doi = {10.1080/00102208108946970}, URL = {https://doi.org/10.1080/00102208108946970}, year = {1981}} +@article{zhong2023, + author = {Hongtao Zhong and Xingqian Mao and Ning Liu and Ziyu Wang and Timothy Ombrello and Yiguang Ju}, + doi = {10.1016/j.combustflame.2023.112948}, + issn = {15562921}, + pages = {112948}, + journal = {Combust. Flame}, + keywords = {Ammonia,Laser diagnostics,N2O/NOx chemistry,Nitrogen cycle,Plasmas,Sustainable energy}, + month = {10}, + publisher = {Elsevier Inc.}, + title = {{Understanding non-equilibrium N$_2$O/NO$_x$ chemistry in plasma-assisted low-temperature NH$_3$ oxidation}}, + volume = {256}, + year = {2023} +} @article{zhu2006, author = {H.~Zhu and R.~J.~Kee}, journal = {Journal of The Electrochemical Society}, diff --git a/doc/sphinx/yaml/reactions.md b/doc/sphinx/yaml/reactions.md index 93c98cff29..44f80ee7b7 100644 --- a/doc/sphinx/yaml/reactions.md +++ b/doc/sphinx/yaml/reactions.md @@ -18,6 +18,7 @@ The fields common to all `reaction` entries are: - [`Blowers-Masel`](sec-yaml-Blowers-Masel) - [`two-temperature-plasma`](sec-yaml-two-temperature-plasma) - [`electron-collision-plasma`](sec-yaml-electron-collision-plasma) + - [`vibrational-relaxation`](sec-yaml-vibrational-relaxation) - [`electron-collisions`](sec-yaml-electron-collisions) - [`falloff`](sec-yaml-falloff) - [`chemically-activated`](sec-yaml-chemically-activated) @@ -118,20 +119,262 @@ a mapping with fields: : The pre-exponential factor `b` -: The temperature exponent, which is applied to the electron temperature +: The electron temperature exponent + +`b-gas` +: The gas temperature exponent `Ea-gas` : The activation energy term $E_{a,g}$ that is related to the gas temperature `Ea-electron` : The activation energy term $E_{a,e}$ that is related to the electron temperature +or a corresponding four-element list. -or a corresponding four-element list. The following are equivalent: +`T-inv` +: Optional gas temperature scale for the term $\exp(-T/T_\mathrm{inv})$. +If omitted or set to 0, this term is not included. + +The following are equivalent: ```yaml -{A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol} +{A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol, b-gas: 0, T-inv: 0} [17283, -3.1, -5820 J/mol, 1081 J/mol] ``` +(sec-vibrational-relaxation-constant-rate)= + +### Vibrational relaxation: constant + +The `constant` vibrational relaxation rate expression is + +$$ +k_f = A +$$ + +where $A$ is a temperature-independent rate coefficient. + +The rate parameters are specified as a mapping with fields: + +`A` +: The pre-exponential factor. The units are the standard Cantera rate coefficient +units for the reaction order and are converted using the unit system. + +No other rate parameters are allowed for this model. In particular, `b`, `n`, `B`, +`C`, `D`, `m`, `E`, `z`, and `scaling` are not permitted. + +This model describes relaxation of one vibrationally excited reactant to its ground +state by collision with an unchanged collider: + +```yaml +- equation: N2(v) + N2 => N2 + N2 + type: vibrational-relaxation + vibration-model: constant + rate-constant: {A: 1.0e2 cm^3/mol/s} +``` + +(sec-vibrational-relaxation-multi-state-resolved-rate)= + +### Vibrational relaxation: multi-state-resolved + +The `multi-state-resolved` vibrational relaxation rate expression is + +$$ +k_f = +\mathrm{scaling} , A +\exp \left( +b \ln T ++ B ++ C T^{-1/3} ++ D T^{-2/3} +\right) +$$ + +where $T$ is the gas temperature in K. + +The rate parameters are specified as a mapping with fields: + +`A` +: The pre-exponential factor. The units are the standard Cantera rate coefficient +units for the reaction order and are converted using the unit system. + +`b` +: The dimensionless temperature exponent. Defaults to 0.0. + +`B` +: Dimensionless constant term in the exponential. Defaults to 0.0. + +`C` +: Coefficient multiplying $T^{-1/3}$. This value is interpreted assuming that +$T$ is in K. Defaults to 0.0. + +`D` +: Coefficient multiplying $T^{-2/3}$. This value is interpreted assuming that +$T$ is in K. Defaults to 0.0. + +`scaling` +: Dimensionless multiplicative scaling factor applied to the rate coefficient. +Defaults to 1.0. + +The parameters `n`, `m`, `E`, and `z` are not permitted for this model. + +This model is intended for detailed V-T and V-V relaxation reactions involving +vibrationally resolved species. All vibrational species in the reaction must belong +to the same vibrational family. + +Example: + +```yaml +- equation: N2(v2) + O => N2(v1) + O + type: vibrational-relaxation + vibration-model: multi-state-resolved + rate-constant: + A: 6.02e+23 cm^3/mol/s + b: 1.0 + B: -34.03 + C: 33.11 + D: 0.0 + scaling: 2.0 +``` + +(sec-vibrational-relaxation-castela-rate)= + +### Vibrational relaxation: castela + +The `castela` vibrational relaxation rate expression is based on the relaxation time + +$$ +\tau_k = +\frac{p_0}{p_k} +\exp \left[ +a_k \left(T^{-1/3} - b_k \right) - 18.42 +\right] +$$ + +and is converted internally to the equivalent bimolecular rate coefficient + +$$ +k_k(T) = +\frac{R T}{p_0} +\exp \left[ +18.42 + a_k b_k - a_k T^{-1/3} +\right] +$$ + +where $T$ is the gas temperature in K, $R$ is the gas constant, and $p_0$ is the +reference pressure. + +The rate parameters are specified as a mapping with fields: + +`a` +: Castela coefficient $a_k$. This value is interpreted assuming that $T$ is in K. + +`b` +: Castela coefficient $b_k$. This value is interpreted assuming that $T$ is in K. + +`reference-pressure` +: Reference pressure $p_0$. Defaults to 1 atm. + +The parameter `A` is not specified by the user for this model. The internal +pre-exponential factor is computed as $R / p_0$. The parameters `A`, `n`, `K`, `B`, +`C`, `D`, `m`, `E`, `z`, and `scaling` are not permitted. + +The `castela` model is only valid for N2 vibrational relaxation and only supports +the colliders `N2`, `O2`, and `O`. It describes relaxation to the ground state, so +the reaction must have the form `N2(v) + M => N2 + M`, where `M` is one of the +supported colliders. + +This model is intended to implement the mean vibrational energy equation model by +the means of a fictitious species N2(v) lumping together all N2 vibrational states. + +Example: + +```yaml +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: castela + rate-constant: + a: 72.4 + b: 0.015 + reference-pressure: 1 atm +``` + +(sec-vibrational-relaxation-starikovskiy-rate)= + +### Vibrational relaxation: starikovskiy + +The `starikovskiy` vibrational relaxation rate expression is + +$$ +k_f = +A T^n +\exp \left( +K ++ B T^{-1/3} ++ C T^{-m} ++ D T^{-z} +\right) +$$ + +where $T$ is the gas temperature in K. + +The rate parameters are specified as a mapping with fields: + +`A` +: The pre-exponential factor. The units are the standard Cantera rate coefficient +units for the reaction order and are converted using the unit system. + +`n` +: The dimensionless temperature exponent. Defaults to 0.0. + +`K` +: Dimensionless constant term in the exponential. Defaults to 0.0. + +`B` +: Signed coefficient multiplying $T^{-1/3}$. This value is interpreted assuming +that $T$ is in K. Defaults to 0.0. + +`C` +: Signed coefficient multiplying $T^{-m}$. This value is interpreted assuming +that $T$ is in K. Defaults to 0.0. + +`m` +: Positive exponent used by the `C` term. Defaults to 1.0. + +`D` +: Signed coefficient multiplying $T^{-z}$. This value is interpreted assuming +that $T$ is in K. Defaults to 0.0. + +`z` +: Positive exponent used by the `D` term. Defaults to 1.0. + +The parameters `b` and `scaling` are not permitted for this model. + +The coefficients `B`, `C`, and `D` are signed values. For example, a negative +coefficient multiplying $T^{-1/3}$ should be written directly as `B: -...` in +the YAML input. + +This model describes relaxation of one vibrationally excited reactant to its ground +state by collision with an unchanged collider: + +```yaml +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: starikovskiy + rate-constant: + A: 6.02214076e+23 cm^3/mol/s + n: 1.0 + K: -34.03 + B: -33.11 + C: 0.0 + m: 1.0 + D: 0.0 + z: 1.0 +``` + +This model is intended to extend the mean vibrational energy equation model from castela +by allowing collisions with more potential colliders. It should be employend alongside a +fictitious species X(v) lumping together all X molecule vibrational states for each molecule +that the user wishes to descibe vibrationally. (sec-yaml-efficiencies)= ## Efficiencies @@ -282,6 +525,51 @@ Example: rate-constant: {A: 17283, b: -3.1, Ea-gas: -5820 J/mol, Ea-electron: 1081 J/mol} ``` +(sec-yaml-vibrational-relaxation)= +### `vibrational-relaxation` +In plasma phases, vibrationally excited species play an important role in +the process of slow gas heating. Their vibrational-translational relaxation +rate is therefore important to model. This class offers four options to +model VT-relaxation rates: +- [`constant`](sec-vibrational-relaxation-constant-rate) +- [`multi-state-resolved`](sec-vibrational-relaxation-multi-state-resolved-rate) +- [`castela`](sec-vibrational-relaxation-castela-rate) +- [`starikovskiy`](sec-vibrational-relaxation-starikovskiy-rate) + +Example: + +```yaml +- equation: N2(v2) + O => N2(v1) + O + type: vibrational-relaxation + vibration-model: multi-state-resolved + rate-constant: + A: 6.02e+23 + b: 1.0 + B: -34.03 + C: 33.11 + D: 0.0 + scaling: 2.0 + +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: constant + rate-constant: {A: 1e2} + +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: castela + rate-constant: + a: 72.4 + b: 0.015 + reference-pressure: 1 atm + +- equation: N2(v) + O => N2 + O + type: vibrational-relaxation + vibration-model: starikovskiy + rate-constant: {A: 6.0221407600e+23, n: 1.0, K: -34.03, B: -33.11, C: 0.0, m: 1.0, D: 0.0, z: 1.0} + +``` + (sec-yaml-electron-collision-plasma)= ### `electron-collision-plasma` diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h index 7f17e9de72..2449af1efe 100644 --- a/include/cantera/kinetics/TwoTempPlasmaRate.h +++ b/include/cantera/kinetics/TwoTempPlasmaRate.h @@ -12,7 +12,7 @@ namespace Cantera { -//! Data container holding shared data specific to TwoTempPlasmaRate +//! Data container holding shared data specific to TwoTempPlasmaRate. /** * The data container `TwoTempPlasmaData` holds precalculated data common to * all `TwoTempPlasmaRate` objects. @@ -38,7 +38,6 @@ struct TwoTempPlasmaData : public ReactionData double recipTe = 1.0; //!< inverse of electron temperature }; - //! Two temperature plasma reaction rate type depends on both //! gas temperature and electron temperature. /*! @@ -48,12 +47,19 @@ struct TwoTempPlasmaData : public ReactionData * activation energy for electron is included. * * @f[ - * k_f = A T_e^b \exp (-E_{a,g}/RT) \exp (E_{a,e} (T_e - T)/(R T T_e)) + * k_f = A T^{b_g} T_e^b + * \exp\left(-\frac{E_{a,g}}{RT}\right) + * \exp\left(\frac{E_{a,e}(T_e - T)}{R T T_e}\right) + * \exp\left(-\frac{T}{T_\mathrm{inv}}\right) * @f] * * where @f$ T_e @f$ is the electron temperature, @f$ E_{a,g} @f$ is the activation * energy for gas, and @f$ E_{a,e} @f$ is the activation energy for electron, see * Kossyi, et al. @cite kossyi1992. + * The optional gas temperature exponent b_g defaults to zero, which strictly corresponds to @cite kossyi1992. + * If b_g is non-zero, a generalisation is used. The optional temperature scale + * @f$ T_\mathrm{inv} @f$ adds the term @f$ \exp(-T/T_\mathrm{inv}) @f$ @cite capitelli2013. + * If @f$ T_\mathrm{inv} = 0 @f$, this term is omitted. * * @ingroup arrheniusGroup */ @@ -66,49 +72,102 @@ class TwoTempPlasmaRate : public ArrheniusBase /*! * @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] + * @param b Electron temperature exponent (non-dimensional). + * @param Ea Activation energy in energy units [J/kmol]. + * @param EE Activation electron energy in energy units [J/kmol]. + * @param bg Gas temperature exponent (non-dimensional). If not specified, defaults to 0. + * @since New in %Cantera 4.0 + */ + TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg); + + //! 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 Electron temperature exponent (non-dimensional). + * @param Ea Activation energy in energy units [J/kmol]. Defaults to 0. + * @param EE Activation electron energy in energy units [J/kmol]. Defaults to 0. */ TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0); + //! Constructor in the cases where T_inv is necessary. + /*! + * @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 Electron temperature exponent (non-dimensional). + * @param Ea Activation energy in energy units [J/kmol]. + * @param EE Activation electron energy in energy units [J/kmol]. Defaults to 0. + * @param bg Gas temperature exponent (non-dimensional). Defaults to 0. + * @param Tinv Temperature scale for the term @f$ \exp(-T/T_\mathrm{inv}) @f$ [K]. + * If zero, this term is omitted. Defaults to 0. + * @since New in %Cantera 4.0 + */ + TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg, + double Tinv); + + //! Constructor based on an AnyMap object instead of all parameters directly. TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={}); + //! Creates a new two-temperature-plasma reaction. unique_ptr newMultiRate() const override { return make_unique>(); } + //! Returns the reaction type. const string type() const override { return "two-temperature-plasma"; } + //! Check that the reaction does not have the 'reversible: true' attribute. void setContext(const Reaction& rxn, const Kinetics& kin) override; - //! Evaluate reaction rate + //! Evaluates reaction rate. /*! - * @param shared_data data shared by all reactions of a given type + * @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); + // m_E4_R is the electron activation energy in temperature units. + double logRate = m_bg * shared_data.logT + + 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; + + if (m_Tinv != 0.0) { + logRate += -shared_data.temperature / m_Tinv; + } + + return m_A * std::exp(logRate); } //! Evaluate derivative of reaction rate with respect to temperature - //! divided by reaction rate + //! 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 + * @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] + //! Return the electron activation energy *Ea* [J/kmol]. double activationElectronEnergy() const { return m_E4_R * GasConstant; } + +protected: + //! Get parameters. + void getParameters(AnyMap& node) const override; + + //! Gas temperature exponent. + //! @since New in %Cantera 4.0 + double m_bg = 0.0; + + //! Temperature scale for the optional term @f$ \exp(-T/T_\mathrm{inv}) @f$. + /*! + * A value of zero disables this term. + * @since New in %Cantera 4.0 + */ + double m_Tinv = 0.0; }; } diff --git a/include/cantera/kinetics/VibrationalRelaxationRate.h b/include/cantera/kinetics/VibrationalRelaxationRate.h new file mode 100644 index 0000000000..5065ac72ed --- /dev/null +++ b/include/cantera/kinetics/VibrationalRelaxationRate.h @@ -0,0 +1,311 @@ +//! @file VibrationalRelaxationRate.h +//! Header for vibrational relaxation reaction rates in plasma kinetics. +//! @since New in %Cantera 4.0 +//! +//! This file is part of Cantera. See License.txt in the top-level directory or +//! at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_VIBRATIONALRELAXATIONRATE_H +#define CT_VIBRATIONALRELAXATIONRATE_H + +#include "Arrhenius.h" + +#include + +namespace Cantera +{ +//! Shared temperature data for vibrational relaxation rates. +struct DetailedVibData : public ReactionData +{ + //! Update cached temperature-dependent data. + /** + * @param phase Thermodynamic phase used to retrieve the gas temperature. + * @param kin Kinetics object. Not used here, but required by the + * ReactionData interface. + * @return `true` if the temperature has changed and rates need to be + * recomputed; `false` otherwise. + */ + bool update(const ThermoPhase& phase, const Kinetics& kin) override; + + using ReactionData::update; +}; + +//! Vibrational relaxation reaction rate class definition. +/** + * This class provides a common implementation for several vibrational + * relaxation models: + * - `constant` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + * Internally, all models are mapped to the following generic expression: + * @f[ + * k_f = + * scaling \, A \, + * \exp \left( + * b \ln T + * + B + * + C T^{-1/3} + * + D T^{-m} + * + E T^{-z} + * \right) + * @f] + * where `T` is the gas temperature in K. + * + * The YAML reaction type is: + * @code{.yaml} + * type: vibrational-relaxation + * @endcode + * The selected physical model is specified separately using: + * @code{.yaml} + * type: vibrational-relaxation + * vibration-model: multi-state-resolved + * @endcode + * Accepted values for `vibration-model` are `constant`, + * `multi-state-resolved`, `starikovskiy`, and `castela`. + * + * The `constant` model relaxes the vibrational species with a constant rate + * coefficient. It could just as well be an Arrhenius rate, but the constant + * model is provided for convenience and to avoid confusion with conventional + * Arrhenius rates for the YAML file user. + * + * The `multi-state-resolved` model fully resolves vibrational relaxation by + * taking into account all vibrational species in the phase, for example + * `N2(v=1-8)`, and solves for their V-T and V-V relaxation. The scaling of the + * rates is based on the SSH theory detailed in Chapter 7 of + * @cite capitelli2013. The simplified SSH theory implemented here is based on + * the harmonic oscillator approximation and can be found in equations 18 and + * 19 of @cite guerra2019. The @f$ k_{10} @f$ rates are taken from + * @cite zhong2023, @cite capitelli2013, and @cite starikovskiy2013. + * + * The `castela` model is meant to be used only for N2 vibrational relaxation, + * by collisions with N2, O2, and O exclusively. It implements the mean + * vibrational energy relaxation model using a fictitious Cantera species and + * is based on @cite castela2016. + * + * The `starikovskiy` model is an extension of the Castela model to several + * vibrational species and additional colliders. Many vibrational relaxation + * rates can be found in Table 1 of @cite starikovskiy2013, hence the model + * name. The rates for the vibrational relaxation of NH3 can be found in the + * reaction mechanism provided in the supplementary material of + * @cite zhong2023. More rates for the vibrational relaxation of + * @f$ \mathrm{CH_4} @f$ can be found in @cite popov2016. + * + * Unit conventions: + * + * - `A` uses standard Cantera rate coefficient units. Its units depend on the + * reaction order and are converted by `ArrheniusBase`. + * - `b`, `B`, `m`, `z`, and `scaling` are dimensionless. + * - `C` is interpreted as K^(1/3), assuming `T` is in K. + * - `D` is interpreted as K^m, assuming `T` is in K. + * - `E` is interpreted as K^z, assuming `T` is in K. + * + * The coefficients `B`, `C`, `D`, `E`, `m`, `z`, and `scaling` are read as + * raw floating-point values. They are not converted by Cantera's unit system. + * + * @ingroup arrheniusGroup + * @since New in %Cantera 4.0 + */ +class VibrationalRelaxationRate : public ArrheniusBase +{ +public: + //! Default constructor. + VibrationalRelaxationRate(); + + //! Constructor using the internal generic representation. + /** + * @param A Pre-exponential factor. + * @param B Dimensionless constant in the exponential. + * @param C Coefficient multiplying T^(-1/3). + * @param D Coefficient multiplying T^(-m). + * @param b Dimensionless temperature exponent. + * @param scaling Dimensionless scaling factor. + * @param m Temperature exponent used by the D term. + * @param E Coefficient multiplying T^(-z). + * @param z Temperature exponent used by the E term. + * @since New in %Cantera 4.0 + */ + VibrationalRelaxationRate(double A, double B, double C, double D, + double b, double scaling = 1.0, + double m = 2.0 / 3.0, + double E = 0.0, double z = 1.0); + + //! Constructor based on AnyMap content. + //! @since New in %Cantera 4.0 + explicit VibrationalRelaxationRate(const AnyMap& node, + const UnitStack& rate_units = {}); + + //! Set rate parameters from an AnyMap. + //! @since New in %Cantera 4.0 + void setParameters(const AnyMap& node, const UnitStack& rate_units) override; + + //! Get rate parameters for YAML serialization. + //! @since New in %Cantera 4.0 + void getParameters(AnyMap& node) const override; + + //! Create a rate evaluator for this reaction rate type. + //! @since New in %Cantera 4.0 + unique_ptr newMultiRate() const override { + return make_unique>(); + } + + //! String identifying this reaction rate type. + //! @since New in %Cantera 4.0 + const string type() const override { + return "vibrational-relaxation"; + } + + //! Set context of reaction rate evaluation. + /** + * Vibrational relaxation rates are intended for irreversible + * non-equilibrium plasma reactions. Reversible reactions are rejected + * because the reverse rate cannot be obtained from conventional + * thermochemistry for these models. + */ + //! @since New in %Cantera 4.0 + void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate the forward rate coefficient. + //! @since New in %Cantera 4.0 + double evalFromStruct(const DetailedVibData& shared_data) const { + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_scaling * m_A * std::exp( + m_b * shared_data.logT + + m_B + + m_C * invT13 + + m_D * std::pow(invT, m_m) + + m_E * std::pow(invT, m_z) + ); + } + + //! Evaluate the scaled temperature derivative. + /** + * This returns: + * + * @f[ + * \frac{1}{k_f} \frac{d k_f}{dT} + * = + * \frac{d \ln k_f}{dT} + * @f] + * + * For the internal generic expression, this is: + * + * @f[ + * \frac{b}{T} + * - \frac{C}{3} T^{-4/3} + * - m D T^{-m-1} + * - z E T^{-z-1} + * @f] + */ + //! @since New in %Cantera 4.0 + double ddTScaledFromStruct(const DetailedVibData& shared_data) const; + +private: + //! Dimensionless constant in the exponential. + double m_B = 0.0; + + //! Coefficient multiplying T^(-1/3). + double m_C = 0.0; + + //! Coefficient multiplying T^(-m). + double m_D = 0.0; + + //! Dimensionless scaling factor. + double m_scaling = 1.0; + + //! Temperature exponent used by the D term. + double m_m = 2.0 / 3.0; + + //! Coefficient multiplying T^(-z). + double m_E = 0.0; + + //! Temperature exponent used by the E term. + double m_z = 1.0; + + //! Castela coefficient a. + double m_castela_a = 0.0; + + //! Castela coefficient b. + double m_castela_b = 0.0; + + //! Castela reference pressure. + double m_referencePressure = OneAtm; + + //! Selected vibrational relaxation model. + /** + * Accepted values: + * - `constant` + * - `multi-state-resolved` + * - `starikovskiy` + * - `castela` + */ + string m_vibration_model = "multi-state-resolved"; + + //! YAML variable names. + string m_B_str = "B"; + string m_C_str = "C"; + string m_D_str = "D"; + string m_scaling_str = "scaling"; + string m_m_str = "m"; + string m_E_str = "E"; + string m_z_str = "z"; + string m_reference_pressure_str = "reference-pressure"; + + //! Configure the ArrheniusBase part from an already-converted internal A value. + /** + * This is needed for models such as Castela, where the user-facing YAML does + * not contain a standard Arrhenius A coefficient. + * @since New in %Cantera 4.0 + */ + void configureBaseFromInternalA(const AnyMap& node, const UnitStack& rate_units, + double A, double b); + + //! Configure the ArrheniusBase part from a YAML A value and an explicit b. + /** + * This is needed for models such as `constant` and `starikovskiy`, + * where the YAML does not contain the standard Arrhenius pair A / b. + * @since New in %Cantera 4.0 + */ + void configureBaseFromYamlA(const AnyMap& node, const UnitStack& rate_units, + const AnyValue& A, double b); + + //! Sets parameters + void setGenericParameters(double B, double C, double D, double m, + double E, double z, double scaling); + + //! Sub-function of setParameters relative to the 'constant' model + void setConstantParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + //! Sub-function of setParameters relative to the 'multi-state-resolved' model + void setMultiStateParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + //! Sub-function of setParameters relative to the 'starikovskiy' model + void setStarikovskiyParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + //! Sub-function of setParameters relative to the 'castela' model + void setCastelaParameters(const AnyMap& node, const AnyMap& rateMap, + const UnitStack& rate_units); + + void storePreExponentialFactor(AnyMap& target, double A) const; + + //! Sub-function of getParameters relative to the 'constant' model + void getConstantParameters(AnyMap& node, AnyMap& rateNode) const; + + //! Sub-function of getParameters relative to the 'multi-state-resolved' model + void getMultiStateParameters(AnyMap& rateNode) const; + + //! Sub-function of getParameters relative to the 'starikovskiy' model + void getStarikovskiyParameters(AnyMap& rateNode) const; + + //! Sub-function of getParameters relative to the 'castela' model + void getCastelaParameters(AnyMap& node, AnyMap& rateNode) const; +}; + +} +#endif \ No newline at end of file diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 00907cb0a9..f0803b4eba 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -1025,6 +1025,34 @@ class PlasmaPhase: public IdealGasPhase //! Work array mutable std::vector m_work; + + //! A structure to describe species serving as mean vibrational energy reservoirs + struct VibrationalReservoirSpecies { + size_t reservoirIndex = npos; // index of the vibrational reservoir species + size_t baseSpeciesIndex = npos; // index of the corresponding ground state phase species + }; + + //! Vector of species serving as mean vibrational energy reservoirs + std::vector m_vibrationalReservoirSpecies; + + //! A boolean flag to update vibrational reservoir species + bool m_vibrationalReservoirSpeciesNeedUpdate = true; + + //! Threshold fraction of vibrational reservoirs to their ground state + //! above which they may be at risk of hindering phase chemistry. + double m_vibrationalMoleFractionThreshold = 1e-2; + + //! The absolute mole threshold below which a the chemistry is + //! assumed to be safe from alterations from vibrational reservoirs. + double m_vibrationalAbsoluteMoleFractionThreshold = 1e-20; + + //! Detects species serving as vibrational reservoirs and + //! adds them to m_vibrationalReservoirSpecies. + void updateVibrationalReservoirSpecies(); + + //! A function to check that vibrational reservoir species + //! are not at risk to hinder phase chemistry. + void checkVibrationalReservoirMoleFractions(); }; } diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 1019fe4f4e..cb079dc0ee 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -18,6 +18,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" namespace Cantera { @@ -40,6 +41,11 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); + // VibrationalRelaxationRate evaluator + reg("vibrational-relaxation", [](const AnyMap& node, const UnitStack& rate_units) { + return new VibrationalRelaxationRate(node, rate_units); + }); + // ElectronCollisionPlasmaRate evaluator reg("electron-collision-plasma", [](const AnyMap& node, const UnitStack& rate_units) { return new ElectronCollisionPlasmaRate(node, rate_units); diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp index ce3f469496..a53ad3256d 100644 --- a/src/kinetics/TwoTempPlasmaRate.cpp +++ b/src/kinetics/TwoTempPlasmaRate.cpp @@ -57,19 +57,60 @@ TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE) m_Ea_str = "Ea-gas"; m_E4_str = "Ea-electron"; m_E4_R = EE / GasConstant; + m_bg = 0.0; + m_Tinv = 0.0; +} + +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, double bg) + : TwoTempPlasmaRate(A, b, Ea, EE, bg, 0.0) +{ +} + +TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE, + double bg, double Tinv) + : TwoTempPlasmaRate(A, b, Ea, EE) +{ + m_bg = bg; + m_Tinv = Tinv; } TwoTempPlasmaRate::TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units) : TwoTempPlasmaRate() { setParameters(node, rate_units); + m_bg = node.getDouble("b-gas", 0.0); + m_Tinv = node.getDouble("T-inv", 0.0); +} + +void TwoTempPlasmaRate::getParameters(AnyMap& node) const +{ + ArrheniusBase::getParameters(node); + + auto& rateNode = node["rate-constant"].as(); + + if (m_bg != 0.0) { + rateNode["b-gas"] = m_bg; + } + + if (m_Tinv != 0.0) { + rateNode["T-inv"] = m_Tinv; + } + + rateNode.setFlowStyle(); } double TwoTempPlasmaRate::ddTScaledFromStruct(const TwoTempPlasmaData& shared_data) const { warn_user("TwoTempPlasmaRate::ddTScaledFromStruct", "Temperature derivative does not consider changes of electron temperature."); - return (m_Ea_R - m_E4_R) * shared_data.recipT * shared_data.recipT; + double derivative = m_bg * shared_data.recipT + (m_Ea_R - m_E4_R) + * shared_data.recipT * shared_data.recipT; + + if (m_Tinv != 0.0) { + derivative += -1.0 / m_Tinv; + } + + return derivative; } void TwoTempPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) diff --git a/src/kinetics/VibrationalRelaxationRate.cpp b/src/kinetics/VibrationalRelaxationRate.cpp new file mode 100644 index 0000000000..2016b26277 --- /dev/null +++ b/src/kinetics/VibrationalRelaxationRate.cpp @@ -0,0 +1,846 @@ +//! @file VibrationalRelaxationRate.cpp +//! @since New in %Cantera 4.0 + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/kinetics/VibrationalRelaxationRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/base/global.h" + +#include +#include +#include +#include +#include +#include + +namespace Cantera +{ + +namespace +{ + +constexpr double VibTolerance = 1e-12; + +const string WhereSetParameters = "VibrationalRelaxationRate::setParameters"; +const string WhereGetParameters = "VibrationalRelaxationRate::getParameters"; +const string WhereSetContext = "VibrationalRelaxationRate::setContext"; + +const string ModelConstant = "constant"; +const string ModelMultiState = "multi-state-resolved"; +const string ModelStarikovskiy = "starikovskiy"; +const string ModelCastela = "castela"; + +enum class VibModel { + Constant, + MultiStateResolved, + Starikovskiy, + Castela +}; + +VibModel parseVibrationModel(const string& model, const AnyMap& input, + const string& where) +{ + if (model == ModelConstant) { + return VibModel::Constant; + } else if (model == ModelMultiState) { + return VibModel::MultiStateResolved; + } else if (model == ModelStarikovskiy) { + return VibModel::Starikovskiy; + } else if (model == ModelCastela) { + return VibModel::Castela; + } + + throw InputFileError(where, input, + "Unrecognized vibration-model '{}'. Expected 'multi-state-resolved', " + "'starikovskiy', 'castela', or 'constant'.", model); +} + +const AnyMap& getRateConstantMap(const AnyMap& node) +{ + if (!node.hasKey("rate-constant")) { + throw InputFileError(WhereSetParameters, node, + "A vibrational-relaxation reaction requires a 'rate-constant' " + "mapping."); + } + + const auto& rate = node["rate-constant"]; + + if (!rate.is()) { + throw InputFileError(WhereSetParameters, node, + "The 'rate-constant' field must be a mapping."); + } + + return rate.as(); +} + +// helpers to check the correctness of a chosen vibrational model input data. +void requireNoKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (node.hasKey(key)) { + throw InputFileError(where, node, + "Key '{}' is not allowed for vibration-model '{}'.", key, model); + } +} + +void requireKey(const AnyMap& node, const string& key, + const string& model, const string& where) +{ + if (!node.hasKey(key)) { + throw InputFileError(where, node, + "Missing required key '{}' for vibration-model '{}'.", key, model); + } +} + +void requireKeys(const AnyMap& node, const string& model, + const string& where, std::initializer_list keys) +{ + for (const auto& key : keys) { + requireKey(node, key, model, where); + } +} + +void forbidKeys(const AnyMap& node, const string& model, + const string& where, std::initializer_list keys) +{ + for (const auto& key : keys) { + requireNoKey(node, key, model, where); + } +} + +// check whether a species is a vibrational excited species +// or a vibrational reservoir from its name +bool isVibrationalSpecies(const string& name) +{ + // Supported names formats: + // + // N2(v) + // O2(v0) + // O2(v1) + // O2(v12) + // + // Style like O2(v=1) should not be used. + const auto pos = name.find("(v"); + return pos != string::npos && !name.empty() && name.back() == ')'; +} + +// find the name of the ground state of a vibrationally excited species. +string groundStateName(const string& name) +{ + const auto pos = name.find("(v"); + if (pos == string::npos) { + return name; + } + return name.substr(0, pos); +} + +string vibrationalFamilyName(const string& name) +{ + // Collapse O2(v1), O2(v2), O2(v12), and O2(v) into O2(v) as a "family + // name" to identify which species have the same ground state "parent". + return groundStateName(name) + "(v)"; +} + +double compositionSum(const Composition& comp) +{ + double sum = 0.0; + for (const auto& item : comp) { + sum += item.second; + } + return sum; +} + +Composition replaceVibrationalSpeciesByGroundState(const Composition& comp) +{ + Composition out; + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (isVibrationalSpecies(name)) { + out[groundStateName(name)] += value; + } else { + out[name] += value; + } + } + return out; +} + +bool sameComposition(Composition diff, const Composition& b, + double tol = VibTolerance) +{ + for (const auto& [name, value] : b) { + diff[name] -= value; + } + + for (const auto& [name, value] : diff) { + if (std::abs(value) > tol) { + return false; + } + } + + return true; +} + +std::vector vibrationalSpeciesInComposition(const Composition& comp) +{ + std::vector out; + + for (const auto& item : comp) { + const string& name = item.first; + const double value = item.second; + + if (value != 0.0 && isVibrationalSpecies(name)) { + out.push_back(name); + } + } + + return out; +} + +// Registry used to check that a given vibrational family uses exactly one +// relaxation model inside one Kinetics object. +// +// It is allowed for each "vibrational family" +// to have its own vibrational relaxation model: +// +// N2(v) -> constant +// O2(v) -> multi-state-resolved +// NH3(v) -> starikovskiy +// +// But it is forbidden for each reaction of the same "vibrational family" +// to have different relaxation models: +// +// N2(v) + O -> castela +// N2(v) + N2 -> starikovskiy +// +std::map> s_modelByFamily; +std::set s_warnedMixedModels; + +void registerVibrationalModelConsistency(const Kinetics& kin, + const string& family, + const string& model, + const AnyMap& input) +{ + auto& modelByFamily = s_modelByFamily[&kin]; + + const auto existing = modelByFamily.find(family); + if (existing != modelByFamily.end() && existing->second != model) { + throw InputFileError(WhereSetContext, input, + "Inconsistent vibration-model for vibrational family '{}'. " + "This family was already registered with model '{}', but the " + "current reaction uses model '{}'. A given vibrational family " + "must use exactly one relaxation model.", + family, existing->second, model); + } + + modelByFamily[family] = model; + + std::set models; + for (const auto& item : modelByFamily) { + models.insert(item.second); + } + + if (models.size() > 1 && !s_warnedMixedModels.count(&kin)) { + std::ostringstream msg; + msg << "Multiple vibrational relaxation models were detected in the " + << "same kinetics object. This is allowed only if each " + << "vibrational family is internally consistent:\n"; + + for (const auto& item : modelByFamily) { + msg << " - " << item.first << ": " << item.second << "\n"; + } + + warn_user(WhereSetContext, msg.str()); + s_warnedMixedModels.insert(&kin); + } +} + +string inferRelaxingFamily(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + + if (vibReactants.empty()) { + throw InputFileError(WhereSetContext, rxn.input, + "A vibrational-relaxation reaction must contain at least one " + "vibrational reactant, for example O2(v1), N2(v), or NH3(v1)."); + } + + return vibrationalFamilyName(vibReactants.front()); +} + +void validateSimpleRelaxationToGroundState(const Reaction& rxn, + const string& model) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.size() != 1) { + throw InputFileError(WhereSetContext, rxn.input, + "vibration-model '{}' expects exactly one vibrational reactant.", + model); + } + + if (!vibProducts.empty()) { + throw InputFileError(WhereSetContext, rxn.input, + "vibration-model '{}' describes relaxation to the ground state " + "and therefore does not allow vibrational products.", model); + } + + const Composition relaxedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + if (!sameComposition(relaxedReactants, rxn.products)) { + throw InputFileError(WhereSetContext, rxn.input, + "Invalid vibrational relaxation stoichiometry for model '{}'. " + "Expected a reaction equivalent to X(v) + M => X + M, where the " + "collider M is unchanged.", model); + } + + if (std::abs(compositionSum(rxn.reactants) - 2.0) > VibTolerance + || std::abs(compositionSum(rxn.products) - 2.0) > VibTolerance) + { + throw InputFileError(WhereSetContext, rxn.input, + "vibration-model '{}' expects a bimolecular relaxation reaction " + "of the form X(v) + M => X + M.", model); + } +} + +void validateCastelaReaction(const Reaction& rxn) +{ + validateSimpleRelaxationToGroundState(rxn, ModelCastela); + + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const string family = vibrationalFamilyName(vibReactants.front()); + + if (family != "N2(v)") { + throw InputFileError(WhereSetContext, rxn.input, + "The Castela relaxation model is only valid for N2 vibrational " + "relaxation. Found vibrational family '{}'.", family); + } + + // Castela parameters used here are only intended for the following + // colliders: N2, O2, and O. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + string collider = ""; + for (const auto& item : collapsedReactants) { + const string& name = item.first; + const double value = item.second; + + if (name == "N2") { + if (std::abs(value - 2.0) < VibTolerance) { + collider = "N2"; + } + } else if (std::abs(value - 1.0) < VibTolerance) { + collider = name; + } + } + + if (collider != "N2" && collider != "O2" && collider != "O") { + throw InputFileError(WhereSetContext, rxn.input, + "The Castela relaxation model only supports colliders 'N2', " + "'O2', and 'O'. Found collider '{}'.", collider); + } +} + +void validateSameVibrationalFamily(const std::vector& species, + const string& family, + const AnyMap& input, + const string& role) +{ + for (const auto& sp : species) { + const string spFamily = vibrationalFamilyName(sp); + if (spFamily != family) { + throw InputFileError(WhereSetContext, input, + "Invalid detailed vibrational relaxation reaction: all " + "vibrational {} must belong to the same vibrational " + "family. Found '{}' and '{}'.", role, family, spFamily); + } + } +} + +void validateDetailedRelaxationReaction(const Reaction& rxn) +{ + const auto vibReactants = vibrationalSpeciesInComposition(rxn.reactants); + const auto vibProducts = vibrationalSpeciesInComposition(rxn.products); + + if (vibReactants.empty()) { + throw InputFileError(WhereSetContext, rxn.input, + "A detailed vibrational relaxation reaction must contain at least " + "one vibrational reactant."); + } + + // All vibrational species in a detailed relaxation reaction are required + // to belong to the same "vibrational family". + const string family = vibrationalFamilyName(vibReactants.front()); + + validateSameVibrationalFamily(vibReactants, family, rxn.input, "reactants"); + validateSameVibrationalFamily(vibProducts, family, rxn.input, "products"); + + // The reaction must conserve the ground-state composition once all + // vibrational labels are removed. + const Composition collapsedReactants = + replaceVibrationalSpeciesByGroundState(rxn.reactants); + + const Composition collapsedProducts = + replaceVibrationalSpeciesByGroundState(rxn.products); + + if (!sameComposition(collapsedReactants, collapsedProducts)) { + throw InputFileError(WhereSetContext, rxn.input, + "Invalid detailed vibrational relaxation reaction: replacing all " + "vibrational species by their ground-state species does not " + "conserve stoichiometry."); + } + + // The current detailed VV/VT formulation is bimolecular. + if (std::abs(compositionSum(rxn.reactants) - 2.0) > VibTolerance + || std::abs(compositionSum(rxn.products) - 2.0) > VibTolerance) + { + throw InputFileError(WhereSetContext, rxn.input, + "Invalid detailed vibrational relaxation reaction: expected a " + "bimolecular reaction with two reactant molecules and two product " + "molecules."); + } +} + +} // end of namespace where all the helpers and safety functions are defined. + +bool DetailedVibData::update(const ThermoPhase& phase, const Kinetics& kin) +{ + const double T = phase.temperature(); + + if (T == temperature) { + return false; + } + + ReactionData::update(T); + return true; +} + +// Default constructor +VibrationalRelaxationRate::VibrationalRelaxationRate() +{ +} + +// Constructor +VibrationalRelaxationRate::VibrationalRelaxationRate( + double A, double B, double C, double D, + double b, double scaling, double m, double E, double z) + : ArrheniusBase(A, b, 0.0) + , m_B(B) + , m_C(C) + , m_D(D) + , m_scaling(scaling) + , m_m(m) + , m_E(E) + , m_z(z) +{ +} + +// AnyMap constructor +VibrationalRelaxationRate::VibrationalRelaxationRate( + const AnyMap& node, const UnitStack& rate_units) + : VibrationalRelaxationRate() +{ + setParameters(node, rate_units); +} + +void VibrationalRelaxationRate::configureBaseFromInternalA( + const AnyMap& node, const UnitStack& rate_units, double A, double b) +{ + // Store the original input and configure reaction-rate units. + // + // We intentionally do not call ArrheniusBase::setParameters here because + // some vibration models do not expose a standard YAML rate-constant with + // both A and b. Castela is the main example. + // This is also done for future class extension compatibility. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + m_A = A; + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // ArrheniusBase::validate checks this flag. + m_valid = true; +} + +void VibrationalRelaxationRate::configureBaseFromYamlA( + const AnyMap& node, const UnitStack& rate_units, + const AnyValue& A, double b) +{ + // Store the original input and configure reaction-rate units first, so + // conversionUnits() is available for A. + ReactionRate::setParameters(node, rate_units); + setRateUnits(rate_units); + + m_negativeA_ok = node.getBool("negative-A", false); + + // Convert the user-facing YAML A value with Cantera's standard + // rate-coefficient unit conversion. + m_A = node.units().convertRateCoeff(A, conversionUnits()); + m_b = b; + + if (m_A != 0.0) { + m_logA = std::log(std::abs(m_A)); + } else { + m_logA = NAN; + } + + // ArrheniusBase::validate checks this flag. + m_valid = true; +} + +void VibrationalRelaxationRate::setParameters(const AnyMap& node, + const UnitStack& rate_units) +{ + // The reaction rate type is always: + // + // type: vibrational-relaxation + // + // The physical model is selected separately by: + // + // vibration-model: constant + // vibration-model: multi-state-resolved + // vibration-model: starikovskiy + // vibration-model: castela + // + // The default model is multi-state-resolved. + + m_vibration_model = node.getString("vibration-model", ModelMultiState); + const auto& rateMap = getRateConstantMap(node); + + switch (parseVibrationModel(m_vibration_model, node, WhereSetParameters)) { + case VibModel::Constant: + setConstantParameters(node, rateMap, rate_units); + break; + case VibModel::MultiStateResolved: + setMultiStateParameters(node, rateMap, rate_units); + break; + case VibModel::Starikovskiy: + setStarikovskiyParameters(node, rateMap, rate_units); + break; + case VibModel::Castela: + setCastelaParameters(node, rateMap, rate_units); + break; + } +} + +void VibrationalRelaxationRate::setConstantParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // Constant model: + // + // k(T) = A + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_b_str, "n", m_B_str, m_C_str, m_D_str, m_m_str, + m_E_str, m_z_str, m_scaling_str}); + + configureBaseFromYamlA(node, rate_units, rateMap[m_A_str], 0.0); + setGenericParameters(0.0, 0.0, 0.0, 2.0 / 3.0, 0.0, 1.0, 1.0); +} + +void VibrationalRelaxationRate::setMultiStateParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // Detailed VV/VT model: + // + // k(T) = scaling * A * exp( + // b * log(T) + // + B + // + C * T^(-1/3) + // + D * T^(-2/3) + // ) + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {"n", m_m_str, m_E_str, m_z_str}); + + ArrheniusBase::setParameters(node, rate_units); + + setGenericParameters( + rateMap.getDouble(m_B_str, 0.0), + rateMap.getDouble(m_C_str, 0.0), + rateMap.getDouble(m_D_str, 0.0), + 2.0 / 3.0, + 0.0, + 1.0, + rateMap.getDouble(m_scaling_str, 1.0)); +} + +void VibrationalRelaxationRate::setStarikovskiyParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // User-facing formula: + // + // k(T) = A * T^n * exp( + // K + // + B * T^(-1/3) + // + C * T^(-m) + // + D * T^(-z) + // ) + // + // B, C, and D are signed coefficients read directly from YAML. + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {m_A_str}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_b_str, m_scaling_str}); + + const double m = rateMap.getDouble("m", 1.0); + const double z = rateMap.getDouble("z", 1.0); + + if (m <= 0.0 || z <= 0.0) { + throw InputFileError(WhereSetParameters, node, + "The Starikovskiy exponents 'm' and 'z' must be positive."); + } + + configureBaseFromYamlA( + node, rate_units, rateMap[m_A_str], rateMap.getDouble("n", 0.0)); + + setGenericParameters( + rateMap.getDouble("K", 0.0), + rateMap.getDouble("B", 0.0), + rateMap.getDouble("C", 0.0), + m, + rateMap.getDouble("D", 0.0), + z, + 1.0); +} + +void VibrationalRelaxationRate::setCastelaParameters( + const AnyMap& node, const AnyMap& rateMap, const UnitStack& rate_units) +{ + // Castela model: + // + // Original relaxation time: + // + // tau_k = p0 / p_k + // * exp[a_k * (T^(-1/3) - b_k) - 18.42] + // + // Equivalent bimolecular rate coefficient: + // + // k_k(T) = R T / p0 + // * exp[18.42 + a_k b_k - a_k T^(-1/3)] + // + // Internal mapping: + // + // A = R / p0 + // b = 1 + // B = 18.42 + a_k b_k + // C = -a_k + // D = 0 + // E = 0 + // scaling = 1 + requireKeys(rateMap, m_vibration_model, WhereSetParameters, {"a", "b"}); + + forbidKeys(rateMap, m_vibration_model, WhereSetParameters, + {m_A_str, "n", "K", m_B_str, m_C_str, m_D_str, m_m_str, + m_E_str, m_z_str, m_scaling_str}); + + m_castela_a = rateMap["a"].asDouble(); + m_castela_b = rateMap["b"].asDouble(); + + if (rateMap.hasKey(m_reference_pressure_str)) { + m_referencePressure = rateMap.convert(m_reference_pressure_str, "Pa"); + } else { + m_referencePressure = OneAtm; + } + + if (m_referencePressure <= 0.0) { + throw InputFileError(WhereSetParameters, node, + "Castela reference-pressure must be positive."); + } + + configureBaseFromInternalA( + node, rate_units, GasConstant / m_referencePressure, 1.0); + + setGenericParameters( + 18.42 + m_castela_a * m_castela_b, + -m_castela_a, + 0.0, + 2.0 / 3.0, + 0.0, + 1.0, + 1.0); +} + +void VibrationalRelaxationRate::getParameters(AnyMap& node) const +{ + if (!valid()) { + return; + } + + if (allowNegativePreExponentialFactor()) { + node["negative-A"] = true; + } + + node["vibration-model"] = m_vibration_model; + + AnyMap rateNode; + + switch (parseVibrationModel(m_vibration_model, node, WhereGetParameters)) { + case VibModel::Constant: + getConstantParameters(node, rateNode); + break; + case VibModel::MultiStateResolved: + getMultiStateParameters(rateNode); + break; + case VibModel::Starikovskiy: + getStarikovskiyParameters(rateNode); + break; + case VibModel::Castela: + getCastelaParameters(node, rateNode); + break; + } + + rateNode.setFlowStyle(); + node["rate-constant"] = std::move(rateNode); +} + +void VibrationalRelaxationRate::storePreExponentialFactor( + AnyMap& target, double A) const +{ + if (conversionUnits().factor() != 0.0) { + target[m_A_str].setQuantity(A, conversionUnits()); + } else { + target[m_A_str] = A; + target["__unconvertible__"] = true; + } +} + +void VibrationalRelaxationRate::getConstantParameters( + AnyMap& node, AnyMap& rateNode) const +{ + if (std::abs(m_b) > VibTolerance + || std::abs(m_B) > VibTolerance + || std::abs(m_C) > VibTolerance + || std::abs(m_D) > VibTolerance + || std::abs(m_E) > VibTolerance) + { + throw InputFileError(WhereGetParameters, node, + "Cannot serialize this rate as 'constant': the internal " + "parameters contain temperature-dependent terms."); + } + + storePreExponentialFactor(rateNode, m_scaling * m_A); +} + +void VibrationalRelaxationRate::getMultiStateParameters( + AnyMap& rateNode) const +{ + storePreExponentialFactor(rateNode, m_A); + + rateNode[m_b_str] = m_b; + rateNode[m_B_str] = m_B; + rateNode[m_C_str] = m_C; + rateNode[m_D_str] = m_D; + rateNode[m_scaling_str] = m_scaling; +} + +void VibrationalRelaxationRate::getStarikovskiyParameters( + AnyMap& rateNode) const +{ + storePreExponentialFactor(rateNode, m_A); + + rateNode["n"] = m_b; + rateNode["K"] = m_B; + rateNode["B"] = m_C; + rateNode["C"] = m_D; + rateNode["m"] = m_m; + rateNode["D"] = m_E; + rateNode["z"] = m_z; +} + +void VibrationalRelaxationRate::getCastelaParameters( + AnyMap& node, AnyMap& rateNode) const +{ + if (std::abs(m_b - 1.0) > VibTolerance + || std::abs(m_D) > VibTolerance + || std::abs(m_E) > VibTolerance + || std::abs(m_scaling - 1.0) > VibTolerance) + { + throw InputFileError(WhereGetParameters, node, + "Cannot serialize this rate as 'castela': the internal " + "parameters are not consistent with the Castela form. " + "Expected b = 1, D = 0, E = 0, and scaling = 1."); + } + + if (m_referencePressure <= 0.0) { + throw InputFileError(WhereGetParameters, node, + "Cannot serialize this rate as 'castela': " + "reference-pressure must be positive."); + } + + rateNode["a"] = m_castela_a; + rateNode["b"] = m_castela_b; + rateNode[m_reference_pressure_str].setQuantity(m_referencePressure, "Pa"); +} + +double VibrationalRelaxationRate::ddTScaledFromStruct( + const DetailedVibData& shared_data) const +{ + const double invT = shared_data.recipT; + const double invT13 = std::cbrt(invT); + + return m_b * invT + - (m_C / 3.0) * invT13 * invT + - m_m * m_D * std::pow(invT, m_m) * invT + - m_z * m_E * std::pow(invT, m_z) * invT; +} + +void VibrationalRelaxationRate::setContext(const Reaction& rxn, const Kinetics& kin) +{ + if (rxn.reversible) { + throw InputFileError(WhereSetContext, rxn.input, + "Vibrational relaxation rates do not support reversible " + "reactions."); + } + + const string family = inferRelaxingFamily(rxn); + + switch (parseVibrationModel(m_vibration_model, rxn.input, WhereSetContext)) { + case VibModel::Constant: + validateSimpleRelaxationToGroundState(rxn, ModelConstant); + break; + case VibModel::Castela: + validateCastelaReaction(rxn); + break; + case VibModel::Starikovskiy: + validateSimpleRelaxationToGroundState(rxn, ModelStarikovskiy); + break; + case VibModel::MultiStateResolved: + validateDetailedRelaxationReaction(rxn); + break; + } + + registerVibrationalModelConsistency( + kin, family, m_vibration_model, rxn.input); +} + +void VibrationalRelaxationRate::setGenericParameters( + double B, double C, double D, double m, double E, double z, double scaling) +{ + m_B = B; + m_C = C; + m_D = D; + m_m = m; + m_E = E; + m_z = z; + m_scaling = scaling; +} + +} // namespace Cantera \ No newline at end of file diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 48eeb01bfe..11fd3f7733 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -19,6 +19,25 @@ namespace Cantera { namespace { const double gamma = sqrt(2 * ElectronCharge / ElectronMass); + + //! A function to check whether a phase species is actually a vibrational reservoir species + //! If the species is actually a reservoir species, the function sets the base name + //! at the given baseName adress. + bool isVibrationalReservoirName(const string& name, string& baseName) + { + const string suffix = "(v)"; + + if (name.size() <= suffix.size()) { + return false; + } + + if (name.compare(name.size() - suffix.size(), suffix.size(), suffix) != 0) { + return false; + } + + baseName = name.substr(0, name.size() - suffix.size()); + return !baseName.empty(); + } } PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) @@ -109,6 +128,13 @@ bool PlasmaPhase::addSpecies(shared_ptr spec) "Only one electron species is allowed.", spec->name); } } + + // Adding species may introduce new `X(v)` reservoirs or their base species. + // Set the flag to true to run a check. + if (added) { + m_vibrationalReservoirSpeciesNeedUpdate = true; + } + return added; } @@ -1006,6 +1032,91 @@ double PlasmaPhase::intrinsicHeating() const double qJ = jouleHeatingPower(); checkFinite(qJ); + // set the check here to be updated at runtime + checkVibrationalReservoirMoleFractions(); + return qJ; } + +void PlasmaPhase::updateVibrationalReservoirSpecies() +{ + m_vibrationalReservoirSpecies.clear(); + + for (size_t k = 0; k < nSpecies(); k++) { + string baseName; + const string& reservoirName = speciesName(k); + + if (!isVibrationalReservoirName(reservoirName, baseName)) { + continue; + } + + size_t kBase = speciesIndex(baseName, false); + + // Check that the base species of the vibrational reservoir actually exists in the phase. + if (kBase == npos) { + warn_user("PlasmaPhase::updateVibrationalReservoirSpecies", + "Species '{}' matches the fictive vibrational-reservoir naming " + "convention, but the corresponding base species '{}' was not found. " + "This species will not be monitored as a vibrational reservoir.", + reservoirName, baseName); + continue; + } + + VibrationalReservoirSpecies reservoir; + reservoir.reservoirIndex = k; + reservoir.baseSpeciesIndex = kBase; + + m_vibrationalReservoirSpecies.push_back(reservoir); + } + // the update being done, the flag can be set back to false. + m_vibrationalReservoirSpeciesNeedUpdate = false; +} + +void PlasmaPhase::checkVibrationalReservoirMoleFractions() +{ + if (m_vibrationalReservoirSpeciesNeedUpdate) { + updateVibrationalReservoirSpecies(); + } + + for (const auto& reservoir : m_vibrationalReservoirSpecies) { + const size_t kReservoir = reservoir.reservoirIndex; + const size_t kBase = reservoir.baseSpeciesIndex; + + const double Xv = moleFraction(kReservoir); // mole fraction of the vibrational reservoir species + const double Xb = moleFraction(kBase); // mole fraction of the base species it is associated to + + const double pool = Xv + Xb; + + // Ignore species pools that are too diluted to meaningfully affect chemistry. + if (pool <= m_vibrationalAbsoluteMoleFractionThreshold) { + continue; + } + + // Check that the fraction of vibrational species is not too high + // with respect to its base species. Should this fraction be too high, there + // is a risk for the phase chemistry to be altered by the reservoir: + // the code raises a warning to the user. + const double reservoirFraction = Xv / pool; + + if (reservoirFraction > m_vibrationalMoleFractionThreshold) { + const string& reservoirName = speciesName(kReservoir); + const string& baseName = speciesName(kBase); + + warn_user("PlasmaPhase::checkVibrationalReservoirMoleFractions", + "Warning: fictive vibrational reservoir species '{}' contains " + "{:.3e} of the total '{}' pool. " + "X({}) = {:.3e}, X({}) = {:.3e}, threshold = {:.3e}. " + "Chemistry involving '{}' may be affected because part of the " + "material is stored in an inert vibrational reservoir.\n", + reservoirName, + reservoirFraction, + baseName, + reservoirName, Xv, + baseName, Xb, + m_vibrationalMoleFractionThreshold, + baseName); + } + } +} + } diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index b4d120099d..332ab9cb11 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -14,6 +14,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/base/Array.h" #include "cantera/base/stringUtils.h" diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index fbe699b388..d03da10128 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -13,6 +13,7 @@ #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" #include "cantera/kinetics/TwoTempPlasmaRate.h" +#include "cantera/kinetics/VibrationalRelaxationRate.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/base/Array.h"