[Test] Add tests for more complex LinearBurkeRate cases

Cantera · Oct 31, 2024 · c8805b7 · c8805b7
1 parent 2f71751
commit c8805b7
Show file tree

Hide file tree

Showing 3 changed files with 162 additions and 0 deletions.
diff --git a/src/kinetics/LinearBurkeRate.cpp b/src/kinetics/LinearBurkeRate.cpp
@@ -292,6 +292,7 @@ double LinearBurkeRate::evalFromStruct(const LinearBurkeData& shared_data)
     // logPeff = shared_data.logP + log(eps_mix) - log(eps_M)
     // but log(eps_M)=0 always
     logPeff = shared_data.logP + log(eps_mix);
+
     // We actually have
     // k_LMR_+=evalPlogRate(shared_data,m_dataObj_M,m_rateObj_M)*eps_M*sigmaX_M/eps_mix
     // k_LMR_+=evalTroeRate(shared_data,m_dataObj_M,m_rateObj_M)*eps_M*sigmaX_M/eps_mix

diff --git a/test/data/linearBurke-test.yaml b/test/data/linearBurke-test.yaml
@@ -35,6 +35,22 @@ phases:
   transport: mixture-averaged
   state: {T: 300.0, P: 1 atm}
 
+- name: linear-Burke-complex-eig0
+  thermo: ideal-gas
+  species:
+  - pdep-test.yaml/species: all
+  kinetics: bulk
+  reactions: [reactions-complex-eig0]
+  state: {T: 300.0, P: 1 atm}
+
+- name: linear-Burke-complex-efficiency
+  thermo: ideal-gas
+  species:
+  - pdep-test.yaml/species: all
+  kinetics: bulk
+  reactions: [reactions-complex-efficiency]
+  state: {T: 300.0, P: 1 atm}
+
 species:
 - name: H
   composition: {H: 1}
@@ -248,6 +264,78 @@ linear-Burke_reactions:
   - name: H2O
     efficiency: {A: 10, b: 0, Ea: 0}
 
+# A more complex linear-Burke definition with multiple colliders that use different
+# parameterizations. Duplicate reactions using the baseline parameterizations are
+# included for rate comparisons
+reactions-complex-efficiency:
+- equation: R1A + R1B (+M) <=> P1 + H (+M) # Reaction 1
+  type: linear-Burke
+  duplicate: true
+  colliders:
+  - name: M
+    type: Chebyshev
+    temperature-range: [300.0, 2000.0]
+    pressure-range: [1e3, 1e7]
+    data:
+    - [8.2883, -1.1397, -0.12059, 0.016034]
+    - [1.9764, 1.0037, 7.2865e-03, -0.030432]
+    - [0.3177, 0.26889, 0.094806, -7.6385e-03]
+    - [-0.031285, -0.039412, 0.044375, 0.014458]
+  - name: P3A
+    efficiency: {A: 3.0, b: 0.0, Ea: 0.0}
+    type: pressure-dependent-Arrhenius
+    rate-constants:
+    - {P: 0.01 atm, A: 1.2124e+16, b: -0.5779, Ea: 1.08727e+04}
+    - {P: 1.0 atm, A: 4.9108e+31, b: -4.8507, Ea: 2.47728e+04}
+    - {P: 10.0 atm, A: 1.2866e+47, b: -9.0246, Ea: 3.97965e+04}
+    - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04}
+  - name: P3B
+    efficiency: {A: 5.0, b: 0.0, Ea: 0.0}
+    type: falloff
+    low-P-rate-constant: {A: 2.3e+18, b: -0.9, Ea: -1700.0}
+    high-P-rate-constant: {A: 7.4e+13, b: -0.37, Ea: 0.0}
+    Troe: {A: 0.7346, T3: 94.0, T1: 1756.0, T2: 5182.0}
+  - name: P5A  # identical rate to M, unity efficiency
+    type: Chebyshev
+    efficiency: {A: 1.0, b: 0.0, Ea: 0.0}
+    temperature-range: [300.0, 2000.0]
+    pressure-range: [1e3, 1e7]
+    data:
+    - [8.2883, -1.1397, -0.12059, 0.016034]
+    - [1.9764, 1.0037, 7.2865e-03, -0.030432]
+    - [0.3177, 0.26889, 0.094806, -7.6385e-03]
+    - [-0.031285, -0.039412, 0.044375, 0.014458]
+  - name: R6
+    efficiency: {A: 7.0, b: 0.0, Ea: 0.0}
+
+- equation: R1A + R1B <=> P1 + H  # M for Reaction 1
+  type: Chebyshev
+  duplicate: true
+  temperature-range: [300.0, 2000.0]
+  pressure-range: [1e3, 1e7]
+  data:
+  - [8.2883, -1.1397, -0.12059, 0.016034]
+  - [1.9764, 1.0037, 7.2865e-03, -0.030432]
+  - [0.3177, 0.26889, 0.094806, -7.6385e-03]
+  - [-0.031285, -0.039412, 0.044375, 0.014458]
+
+- equation: R1A + R1B <=> P1 + H  # collider P3A for Reaction 1
+  type: pressure-dependent-Arrhenius
+  duplicate: true
+  rate-constants:
+  - {P: 0.01 atm, A: 1.2124e+16, b: -0.5779, Ea: 1.08727e+04}
+  - {P: 1.0 atm, A: 4.9108e+31, b: -4.8507, Ea: 2.47728e+04}
+  - {P: 10.0 atm, A: 1.2866e+47, b: -9.0246, Ea: 3.97965e+04}
+  - {P: 100.0 atm, A: 5.9632e+56, b: -11.529, Ea: 5.25996e+04}
+
+- equation: R1A + R1B (+M) <=> P1 + H (+M)  # collider P3B for Reaction 1
+  type: falloff
+  duplicate: true
+  low-P-rate-constant: {A: 2.3e+18, b: -0.9, Ea: -1700.0}
+  high-P-rate-constant: {A: 7.4e+13, b: -0.37, Ea: 0.0}
+  Troe: {A: 0.7346, T3: 94.0, T1: 1756.0, T2: 5182.0}
+
+
 # Invalid reaction definitions to test various input errors
 
 reactions-no-colliders:

diff --git a/test/kinetics/rates.cpp b/test/kinetics/rates.cpp
@@ -176,4 +176,77 @@ TEST(InterfaceReaction, CoverageDependency) {
     EXPECT_NEAR(kf[1], 3.7e20 * exp(-(67.4e6-6e6*0.3)/(GasConstant*T)), 1e-14*kf[1]);
 }
 
+TEST(LinearBurkeRate, RateCombinations)
+{
+    auto sol = newSolution("linearBurke-test.yaml", "linear-Burke-complex-efficiency",
+                           "none");
+    auto& thermo = *sol->thermo();
+    auto& kin = *sol->kinetics();
+    size_t nR = kin.nReactions();
+    ASSERT_EQ(nR, 4);
+    double P0 = 2.0 * OneAtm;
+    double T = 800;
+
+    auto kf = [&](size_t iRxn, double T, double P, const string& X="P1:1.0") {
+        thermo.setState_TPX(T, P, X);
+        vector<double> rates(kin.nReactions());
+        kin.getFwdRateConstants(rates.data());
+        return rates[iRxn];
+    };
+
+    double atol = kf(0, T, P0, "R2: 0.6, R3: 0.4") * 1e-9;
+
+    // For a mixture only containing species treated as M, kf should be the same as
+    // this rate coefficient evaluated at the same P
+    EXPECT_NEAR(kf(0, T, P0, "R2: 0.6, R3: 0.4"), kf(1, T, P0), atol);
+
+    // For a collider that specifies the same rate constant as M, kf should be the same
+    EXPECT_NEAR(kf(0, T, P0, "R2: 0.2, P5A: 0.8"),
+                kf(0, T, P0, "R2: 0.6, R3: 0.4"), atol);
+
+    // For a pure mixture of another collider, the rate should be the same as for that
+    // collider treated as a separate reaction (that is, Peff == P)
+    EXPECT_NEAR(kf(0, T, P0, "P3A: 1.0"), kf(2, T, P0), atol);
+    EXPECT_NEAR(kf(0, T, P0, "P3B: 1.0"), kf(3, T, P0), atol);
+
+    // A binary mixture containing one non-M collider (P3A: X = 0.2, efficiency = 3.0)
+    // and one M collider (R2: X = 0.8, efficiency = 1.0)
+    double X_P3A = 0.2;
+    double X_R2 = 1 - X_P3A;
+    double eps_mix = X_P3A * 3.0 + X_R2 * 1.0;
+    double Peff_P3A = P0 * eps_mix / 3.0;
+    double Peff_R2 = P0 * eps_mix;
+    double Pr_P3A = 3.0 / eps_mix * X_P3A; // reduced pressure
+    double Pr_R2 = 1 / eps_mix * X_R2;
+    EXPECT_NEAR(kf(0, T, P0, "P3A: 0.2, R2: 0.8"),
+                kf(2, T, Peff_P3A) * Pr_P3A + kf(1, T, Peff_R2) * Pr_R2, atol);
+
+    // A binary mixture containing two non-M colliders (P3A: X = 0.3, efficiency = 3.0)
+    // and (P3B: X = 0.7, efficiency = 5)
+    X_P3A = 0.3;
+    double X_P3B = 1 - X_P3A;
+    eps_mix = X_P3A * 3.0 + X_P3B * 5.0;
+    Peff_P3A = P0 * eps_mix / 3.0;
+    double Peff_P3B = P0 * eps_mix / 5.0;
+    Pr_P3A = 3.0 / eps_mix * X_P3A;
+    double Pr_P3B = 5.0 / eps_mix * X_P3B;
+    EXPECT_NEAR(kf(0, T, P0, "P3A: 0.3, P3B: 0.7"),
+                kf(2, T, Peff_P3A) * Pr_P3A + kf(3, T, Peff_P3B) * Pr_P3B, atol);
+
+    // A binary mixture containing two non-M colliders, where the second collider
+    // specifies only an efficiency (P3A: X = 0.4, efficiency = 3.0;
+    // R6: X = 0.6, efficiency = 7). For this collider, Peff == Peff_M, but the reduced
+    // pressure Pr is calculated using the efficiency and mole fraction of that species.
+    X_P3A = 0.4;
+    double X_R6 = 1 - X_P3A;
+    eps_mix = X_P3A * 3.0 + X_R6 * 7.0;
+    Peff_P3A = P0 * eps_mix / 3.0;
+    double Peff_M = P0 * eps_mix;
+    Pr_P3A = 3.0 / eps_mix * X_P3A;
+    double Pr_R6 = 7.0 / eps_mix * X_R6;
+    EXPECT_NEAR(kf(0, T, P0, "P3A: 0.4, R6: 0.6"),
+                kf(2, T, Peff_P3A) * Pr_P3A + kf(1, T, Peff_M) * Pr_R6, atol);
+}
+
+
 }