Incorporating the LMR-R mixture rule into Cantera for the first time

doc/doxygen/cantera.bib

@@ -434,6 +434,17 @@ @online{smith1999
         and W.~C.~Gardiner, Jr. and V.~V.~Lissianski and Z.~Qin},
     url = {http://combustion.berkeley.edu/gri-mech/version30/text30.html},
     year = {1999}}
+@article{singal2024,
+    author = {P.~J.~Singal and J.~Lee and L.~Lei and R.~L.~Speth and M.~P.~Burke},
+    journal = {Proceedings of the Combustion Institute},
+    number = {40},
+    pages = {TBD},
+    title = {Implementation of New Mixture Rules Has a Substantial Impact on
+        Combustion Predictions for H2 and NH3},
+    doi = {TBD},
+    url = {TBD},
+    volume = {TBD},
+    year = {2024}} 
 @article{stewart1989,
     author = {P.~H.~Stewart and C.~W.~Larson and D.~Golden},
     journal = {Combustion and Flame},

doc/example-keywords.txt

@@ -18,9 +18,11 @@ fuel cell
 heat transfer
 ignition delay
 internal combustion engine
+jet-stirred reactor
 kinetics
 Matlab
 mixture
+mixture rule
 multicomponent transport
 multiphase
 non-ideal fluid
@@ -41,6 +43,9 @@ saving output
 sensitivity analysis
 strained flame
 surface chemistry
+shock tube
+species profile
+temperature profile
 thermodynamic cycle
 thermodynamics
 transport

doc/sphinx/reference/kinetics/rate-constants.md

@@ -235,6 +235,61 @@ Chebyshev reactions can be defined in the YAML format using the
 [`Chebyshev`](sec-yaml-Chebyshev) reaction `type`.
 ```
 
+(sec-linear-burke)=
+## Linear Burke Rate Expressions
+
+Linear Burke rate expressions employ the reduced-pressure linear mixture rule (LMR-R). This mixture rule is used to evaluate the rate constants of complex-forming reactions, and is a mole-fraction-weighted sum of the bath gas rate constants (when pure) evaluated at the reduced pressure ($R$) and temperature ($T$) of the mixture.
+
+$$
+k_{\text{LMR-R}}(T,P,\textit{\textbf{X}}) = \sum_{i} k_{i}(T,R_{\text{LMR}})\tilde{X}_{i,\text{LMR}}
+$$
+
+where the reduced pressure, $R$, in its most general form
+
+$$
+R_{\text{LMR}}(T,P,\textit{\textbf{X}}) = \frac{\sum_{i} \Lambda_{0,i}(T)X_i[M]}{\Lambda_{\infty}(T)}
+$$
+
+and the fractional contribution of each component to the reduced pressure, $\tilde{X}_{i}$
+
+$$
+\tilde{X}_{i,\text{LMR}}(T,P,\textit{\textbf{X}})=\frac{\Lambda_{0,i}(T)X_i}{\sum_{j} \Lambda_{0,j}(T)X_j}
+$$
+
+can be cast in terms of the absolute value of the least negative chemically significant eigenvalue of the master equation for the $i^{th}$ collider (when pure) in the low-pressure limit, $\Lambda_{0,i}(T)[M]$, and high-pressure limit, $\Lambda_{\infty}(T)$, and $[M]$ is the total concentration.
+
+Evaluating all rate constants at the reduced pressure ($R$)---instead of the pressure ($P$)---of the mixture takes advantage of the fact that rate constants (and their chemically significant eigenvectors) for different colliders are usually far more similar at the same $R$ than the same $P$. In practice, since rate constants are usually expressed with respect to pressure $P$ (which has units of Pa, Torr, bar, atm, etc.) rather than reduced pressure $R$ (which is dimensionless), one needs to find the effective pressure for the $i^{th}$ collider, $P_{i}^{\text{ eff}}$ (with units of $P$), such that the reduced pressure of pure collider $i$ is equal to the reduced pressure of the mixture, which can be shown to be
+
+$$
+P_{i,\text{LMR}}^{\text{ eff}}(T,P,\textit{\textbf{X}}) = \frac{\sum_{j} \Lambda_{0,j}(T)X_j}{\Lambda_{0,i}(T)}P
+$$
+
+such that an alternate version of the generalized LMR-R equation can be written as
+
+\begin{equation}
+k_{\text{LMR-R}}(T,P,\textit{\textbf{X}}) = \sum_{i} k_{i}(T,P_{i,\text{LMR}}^{\text{ eff}})\tilde{X}_{i,\text{LMR}}
+\label{eq:LMRR_k_P_i_eff}
+\end{equation}
+
+The user can either specify $k_i(T,P)$ and $\Lambda_{0,i}(T)$ for each collider, or specify the third-body efficiency, $\epsilon_{0,i}(T)$, for each non-M collider and assign $\epsilon_{0,\text{M}}(T)=1$.  The pressure-dependent aspect of each i-th collider ($k_i(T,P)$) can be specified in any combination of Troe, PLOG, or Chebyshev formats for colliders for which data are available and $\Lambda_{0,i}(T)$ for each collider (or $\epsilon_{0,i}(T)$ for each non-M collider) can be specified in modified Arrhenius format for colliders for which data are available.
+
+In fact, if $\Lambda_{0,i}(T)$ (or $\epsilon_{0,i}(T)$) for any colliders have available data or can be estimated using typical values (as is typically done in kinetic models for reactions in modified Lindemann expressions) but no data for $k_i(T,P)$ are available, this implementation also allows some colliders to be specified with unique $\Lambda_{0,i}(T)$ (or $\epsilon_{0,i}(T)$) without $k_i(T,P)$, by assuming the same reduced-pressure dependence as M (i.e. $k_{i}(T,R)=k_{M}(T,R)$). Cantera employs the following derived form of the generalized LMR-R equation, where the sum over $n$ is only for the colliders for which unique $k_n(T,P)$ are available. If unique $k_i(T,P)$ data are available for all colliders, then the second term in the above equation effectively disappears.
+
+$$
+k_{\text{LMR-R}}(T,P,\textit{\textbf{X}}) &= \sum_{n} k_{n}(T,P_{n,\text{LMR}}^{\text{ eff}})\tilde{X}_{n,\text{LMR}} + k_{M}(T,P_{M,\text{LMR}}^{\text{ eff}}) \left(1-\sum_{n}\tilde{X}_{n,\text{LMR}}\right)
+$$
+
+If the user has limited or incomplete access to parameter inputs (a likely scenario, given the scarcity of puplished third-body efficiencies and master equation eigenvalues), this computational implementation allowes them much greater flexibility and power to make educated assumptions. Further description of the LMR-R theory and computational method is available in {cite:t}`singal2024`.
+
+```{admonition} YAML Usage
+:class: tip
+Linear Burke rate expressions can be defined in the YAML format using the
+[`linear-burke`](sec-yaml-linear-burke) reaction `type`.
+```
+
+```{versionadded} 3.1
+```
+
 (sec-blowers-masel)=
 ## Blowers-Masel Reactions
 

doc/sphinx/yaml/reactions.rst

@@ -397,6 +397,118 @@ Example::
            [-1.43220e-01,  7.71110e-02,  1.27080e-02, -6.41540e-04]]
 
 
+.. _sec-yaml-linear-burke:
+
+``linear-burke``
+-------------
+
+A complex-forming reaction (one that depends on both P and X) parameterized
+according to the reduced-pressure linear mixture rule as
+:ref:`described here <sec-linear-burke>`.
+
+Additional fields are:
+
+``collider-list``
+    A list of dictionaries, where each entry contains parameters corresponding
+    to individual colliders (species in the bath gas).
+
+``collider``
+    The name of the collider species, which must be entered inside quotations (e.g.,
+    ``"H2O"``). The first collider defined must be ``"M"``, which represents the generic
+    reference collider (often ``Ar`` or ``N2``) that represents all species lacking their
+    own explicit parameterization.
+
+``eps`` or ``eig0``
+    The fractional contribution of each bath gas component (collider) to the reduced
+    pressure. ``eps`` represents the third-body efficiency of the collider relative
+    to that of the reference collider ``"M"`` (``eps: {A:1, b:0, Ea: 0}`` must be provided for
+    ``"M"`` by necessity). ``eig0`` represents the absolute value of the least negative chemically
+    significant eigenvalue of the master equation, evaluated for a collider at its low-pressure
+    limit. All explicitly defined colliders must include either ``eps`` or ``eig0``, but the choice
+    must remain consistent throughout a single reaction (either all colliders are defined with ``eps``,
+    or all are defined with ``eig0``). In both cases, the parameters are entered in Arrhenius format to
+    enable representation of their temperature-dependence.
+
+The pressure-dependent aspect of the rate constant can be parameterized in the user's choice of
+:ref:`Troe <sec-yaml-falloff>`,
+:ref:`pressure-dependent-arrhenius <sec-yaml-pressure-dependent-Arrhenius>`, or
+:ref:`Chebyshev <sec-yaml-Chebyshev>` representations. The same parameters used for a standalone
+Troe, PLOG, or Chebyshev reaction are then inserted directly below ``eps`` or ``eig0`` for a given collider
+(note: Troe cannot be given its own ``efficiencies`` key). At minimum, this treatment must be applied to ``"M"``.
+However, additional colliders may also be given their own Troe, PLOG, or Chebyshev
+parameterization if so desired. Mixing and matching of types within the same reaction is allowed (e.g., a PLOG
+table for ``"M"``, Troe parameters for ``"H2"``, and Chebyshev data for ``"NH3"``).
+
+A mathematical description of this YAML implementation can be found in Eq. 8 of
+:cite:t:`singal2024`.
+
+Examples::
+
+    equation: H + OH <=> H2O
+    type: linear-burke
+    collider-list:
+    - collider: 'M' # N2 is reference collider (Troe format)
+      eps: {A: 1, b: 0, Ea: 0}
+      low-P-rate-constant: {A: 4.530000e+21, b: -1.820309e+00, Ea: 4.987000e+02}
+      high-P-rate-constant: {A: 2.510000e+13, b: 2.329303e-01, Ea: -1.142000e+02}
+      Troe: {A: 9.995044e-01, T3: 1.0e-30, T1: 1.0e+30}
+    - collider: 'AR'
+      eps: {A: 2.20621e-02, b: 4.74036e-01, Ea: -1.13148e+02}
+    - collider: 'H2O'
+      eps: {A: 1.04529e-01, b: 5.50787e-01, Ea: -2.32675e+02}
+
+    equation: H + O2 (+M) <=> HO2 (+M)  # Including "(+M)" is optional
+    type: linear-burke
+    collider-list:
+    - collider: "M" # Argon is reference collider (PLOG format)
+      eps: {A: 1, b: 0, Ea: 0}
+      rate-constants:
+      - {P: 1.316e-02 atm, A: 9.39968e+14, b: -2.14348e+00, Ea: 7.72730e+01}
+      - {P: 1.316e-01 atm, A: 1.07254e+16, b: -2.15999e+00, Ea: 1.30239e+02}
+      - {P: 3.947e-01 atm, A: 3.17830e+16, b: -2.15813e+00, Ea: 1.66994e+02}
+      - {P: 1.000e+00 atm, A: 7.72584e+16, b: -2.15195e+00, Ea: 2.13473e+02}
+      - {P: 3.000e+00 atm, A: 2.11688e+17, b: -2.14062e+00, Ea: 2.79031e+02}
+      - {P: 1.000e+01 atm, A: 6.53093e+17, b: -2.13213e+00, Ea: 3.87493e+02}
+      - {P: 3.000e+01 atm, A: 1.49784e+18, b: -2.10026e+00, Ea: 4.87579e+02}
+      - {P: 1.000e+02 atm, A: 3.82218e+18, b: -2.07057e+00, Ea: 6.65984e+02}
+    - collider: "HE"
+      eps: {A: 3.37601e-01, b: 1.82568e-01, Ea: 3.62408e+01}
+    - collider: "N2"
+      eps: {A: 1.24932e+02, b: -5.93263e-01, Ea: 5.40921e+02}
+    - collider: "H2"
+      eps: {A: 3.13717e+04, b: -1.25419e+00, Ea: 1.12924e+03}
+    - collider: "CO2"
+      eps: {A: 1.62413e+08, b: -2.27622e+00, Ea: 1.97023e+03}
+    - collider: "NH3"
+      eps: {A: 4.97750e+00, b: 1.64855e-01, Ea: -2.80351e+02}
+    - collider: "H2O"
+      eps: {A: 3.69146e+01, b: -7.12902e-02, Ea: 3.19087e+01}
+
+    equation: H2O2 <=> 2 OH
+    type: linear-burke
+    collider-list:
+    - collider: 'M' # Argon is reference collider (Chebyshev format)
+      eps: {A: 1, b: 0, Ea: 0}
+      temperature-range: [200.0, 2000.0]
+      pressure-range: [1.000e-01 atm, 1.000e+02 atm]
+      data:
+      - [-1.5843e+01, 8.7088e-01, -9.4364e-02, -2.8099e-03, -4.4803e-04, 1.5809e-03, -2.5088e-04]
+      - [2.3154e+01, 5.2739e-01, 2.8862e-02, -5.4601e-03, 7.0783e-04, -3.0282e-03, 7.8121e-04]
+      - [-3.8008e-01, 8.6349e-02, 4.0292e-02, -7.2269e-03, 5.7570e-04, 2.7944e-03, -1.4912e-03]
+      - [-1.4800e-01, -7.1798e-03, 2.2052e-02, 6.2269e-03, -5.9801e-03, -8.2205e-06, 1.9243e-03]
+      - [-6.0604e-02, -1.4203e-02, 1.3414e-03, 9.6228e-03, 1.7002e-03, -3.6506e-03, -4.3168e-04]
+      - [-2.4557e-02, -9.7102e-03, -5.8753e-03, 3.0456e-03, 5.8666e-03, 1.5037e-03, -2.0073e-03]
+      - [-1.5400e-02, -5.2427e-03, -6.9148e-03, -5.9440e-03, -1.2183e-03, 2.1694e-03, 1.5925e-03]
+    - collider: 'N2'
+      eps: {A: 1.14813e+00, b: 4.60090e-02, Ea: -2.92413e+00}
+    - collider: 'CO2'
+      eps: {A: 8.98839e+01, b: -4.27974e-01, Ea: 2.41392e+02}
+    - collider: 'H2O2'
+      eps: {A: 6.45295e-01, b: 4.26266e-01, Ea: 4.28932e+01}
+    - collider: 'H2O'
+      eps: {A: 1.36377e+00, b: 3.06592e-01, Ea: 2.10079e+02}
+
+
 .. _sec-yaml-interface-Arrhenius:
 
 ``interface-Arrhenius``

include/cantera/kinetics/LinearBurkeRate.h

@@ -0,0 +1,140 @@
+//! @file LINEARBURKERATE.h
+// 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_LINEARBURKERATE_H
+#define CT_LINEARBURKERATE_H
+#include "cantera/kinetics/Arrhenius.h"
+#include <boost/variant.hpp>
+#include "cantera/kinetics/Falloff.h"
+#include "cantera/kinetics/ChebyshevRate.h"
+#include "cantera/kinetics/PlogRate.h"
+
+namespace Cantera
+{
+
+//! Data container holding shared data specific to LinearBurkeRate
+/**
+ * The data container `LmrData` holds precalculated data common to
+ * all `LinearBurkeRate` objects.
+ */
+struct LmrData : public ReactionData
+{
+    LmrData();
+
+    void update(double T, double P) override
+    {
+        ReactionData::update(T);
+        pressure = P;
+        logP = std::log(P);
+    }
+
+    bool update(const ThermoPhase& phase, const Kinetics& kin) override;
+
+    using ReactionData::update;
+
+    //! Perturb pressure of data container
+    /**
+     * The method is used for the evaluation of numerical derivatives.
+     * @param  deltaP  relative pressure perturbation
+     */
+    void perturbPressure(double deltaP);
+
+    void restore() override;
+
+    virtual void resize(size_t nSpecies, size_t nReactions, size_t nPhases) override
+    {
+        moleFractions.resize(nSpecies, NAN);
+        ready = true;
+    }
+
+    void invalidateCache() override
+    {
+        ReactionData::invalidateCache();
+        pressure = NAN;
+    }
+
+    double pressure = NAN; //!< Pressure
+    double logP = 0.0; //!< Logarithm of pressure
+    bool ready = false; //!< Boolean indicating whether vectors are accessible
+    vector<double> moleFractions;
+    int mf_number;
+
+protected:
+    double m_pressure_buf = -1.0;
+};
+
+
+//! Pressure-dependent and composition-dependent reaction rate calculated
+//! according to the reduced-pressure linear mixture rule (LMR-R) developed
+//! at Columbia University. @cite singal2025 [CITATION NOT YET ACTIVE]
+/*!
+ * [ADD IN THE MATHEMATICAL FORMULA]
+ */
+class LinearBurkeRate final : public ReactionRate
+{
+public:
+    //! Default constructor.
+    LinearBurkeRate() = default;
+
+    //! Constructor from Arrhenius rate expressions at a set of pressures
+    explicit LinearBurkeRate(const std::multimap<double, ArrheniusRate>& rates);
+
+    LinearBurkeRate(const AnyMap& node, const UnitStack& rate_units={});
+
+    unique_ptr<MultiRateBase> newMultiRate() const override {
+        return make_unique<MultiRate<LinearBurkeRate, LmrData>>();
+    }
+
+    //! Identifier of reaction rate type
+    const string type() const override { return "linear-burke"; }
+
+    //! Perform object setup based on AnyMap node information
+    /*!
+     *  @param node  AnyMap containing rate information
+     *  @param rate_units  Unit definitions specific to rate information
+     */
+    void setParameters(const AnyMap& node, const UnitStack& rate_units) override;
+
+    void getParameters(AnyMap& rateNode, const Units& rate_units) const;
+    void getParameters(AnyMap& rateNode) const override {
+        return getParameters(rateNode, Units(0));
+    }
+
+    //! Evaluate overall reaction rate, using Troe/PLOG/Chebyshev to evaluate
+    //! pressure-dependent aspect of the reaction
+    /*!
+     *  @param shared_data  data shared by all reactions of a given type
+     */
+    using RateTypes = boost::variant<PlogRate, TroeRate, ChebyshevRate>;
+    using DataTypes = boost::variant<PlogData, FalloffData, ChebyshevData>;
+    double evalPlogRate(const LmrData& shared_data, DataTypes& dataObj, RateTypes& rateObj);
+    double evalTroeRate(const LmrData& shared_data, DataTypes& dataObj, RateTypes& rateObj);
+    double evalChebyshevRate(const LmrData& shared_data, DataTypes& dataObj, RateTypes& rateObj);
+    double evalFromStruct(const LmrData& shared_data);
+
+    void setContext(const Reaction& rxn, const Kinetics& kin) override;
+
+    void validate(const string& equation, const Kinetics& kin) override;
+
+    vector<size_t> colliderIndices;
+    map<string, AnyMap> colliderInfo;
+    // Third-body collision efficiency objects (eps = eig0_i/eig0_M)
+    // epsObjs1 used for k(T,P,X) and eig0_mix calculation
+    // epsObjs2 used for logPeff calculation
+    vector<ArrheniusRate> epsObjs1;
+    vector<ArrheniusRate> epsObjs2;
+    vector<string> colliderNames;
+    vector<RateTypes> rateObjs;
+    vector<DataTypes> dataObjs;
+    RateTypes rateObj_M;
+    DataTypes dataObj_M;
+    // Third-body collision efficiency object for M (eig0_M/eig0_M = 1)
+    ArrheniusRate epsObj_M;
+    size_t nSpecies;
+    double logPeff_;
+    double eps_mix;
+};
+
+}
+#endif

interfaces/cython/cantera/reaction.pxd

@@ -175,6 +175,10 @@ cdef extern from "cantera/kinetics/PlogRate.h" namespace "Cantera":
         CxxPlogRate(multimap[double, CxxArrheniusRate])
         multimap[double, CxxArrheniusRate] getRates()
 
+cdef extern from "cantera/kinetics/LinearBurkeRate.h" namespace "Cantera":
+    cdef cppclass CxxLinearBurkeRate "Cantera::LinearBurkeRate" (CxxReactionRate):
+        CxxLinearBurkeRate()
+        CxxLinearBurkeRate(CxxAnyMap) except +translate_exception
 
 cdef extern from "cantera/kinetics/ChebyshevRate.h" namespace "Cantera":
     cdef cppclass CxxChebyshevRate "Cantera::ChebyshevRate" (CxxReactionRate):