[Input] Create ElectrochemicalReaction objects from YAML definitions

doc/sphinx/yaml/reactions.rst

@@ -28,8 +28,10 @@ The fields of a ``reaction`` entry are:
     - :ref:`Chebyshev <sec-Chebyshev>`
 
     Reactions on surfaces or edges are automatically treated as
-    :ref:`interface <sec-interface-reaction>` reactions, without the need to
-    specify the ``type``.
+    :ref:`interface <sec-interface-reaction>` reactions, and reactions which
+    involve charge transfer between phases are automatically treated as
+    :ref:`electrochemical <sec-electrochemical-reaction>` reactions, without the
+    need to specify the ``type``.
 
 ``duplicate``
     Boolean indicating whether the reaction is a known duplicate of another
@@ -292,3 +294,27 @@ Example::
     equation: 2 H(s) => H2 + 2 Pt(s)
     rate-constant: {A: 3.7e21 cm^2/mol/s, b: 0, Ea: 67400 J/mol}
     coverage-dependencies: {H(s): {a: 0, m: 0, E: -6000 J/mol}}
+
+
+.. _sec-electrochemical-reaction:
+
+``electrochemical``
+-------------------
+
+Interface reactions involving charge transfer between phases,
+as `described here <https://cantera.org/documentation/dev/doxygen/html/d6/ddd/classCantera_1_1ElectrochemicalReaction.html#details>`_.
+
+Includes the fields of an :ref:`sec-interface-reaction` reaction, plus:
+
+``beta``
+    The symmetry factor for the reaction. Default is 0.5.
+
+``exchange-current-density-formulation``
+    Set to ``true`` if the rate constant parameterizes the exchange current
+    density. Default is ``false``.
+
+Example::
+
+    equation: LiC6 <=> Li+(e) + C6
+    rate-constant: [5.74, 0.0, 0.0]
+    beta: 0.4

src/kinetics/Reaction.cpp

@@ -448,12 +448,11 @@ void setupReaction(Reaction& R, const XML_Node& rxn_node)
     R.reversible = (rev == "true" || rev == "yes");
 }
 
-void setupReaction(Reaction& R, const AnyMap& node)
-{
+void parseReactionEquation(Reaction& R, const AnyValue& equation) {
     // Parse the reaction equation to determine participating species and
     // stoichiometric coefficients
     std::vector<std::string> tokens;
-    tokenizeString(node["equation"].asString(), tokens);
+    tokenizeString(equation.asString(), tokens);
     tokens.push_back("+"); // makes parsing last species not a special case
 
     size_t last_used = npos; // index of last-used token
@@ -478,15 +477,15 @@ void setupReaction(Reaction& R, const AnyMap& node)
                 try {
                     stoich = fpValueCheck(tokens[i-2]);
                 } catch (CanteraError& err) {
-                    throw InputFileError("fpValueCheck", node["equation"],
+                    throw InputFileError("fpValueCheck", equation,
                         err.getMessage());
                 }
             } else {
-                throw InputFileError("setupReaction", node["equation"],
+                throw InputFileError("setupReaction", equation,
                     "Error parsing reaction string '{}'.\n"
                     "Current token: '{}'\nlast_used: '{}'",
-                    node["equation"].asString(),
-                    tokens[i], (last_used == npos) ? "n/a" : tokens[last_used]);
+                    equation.asString(), tokens[i],
+                    (last_used == npos) ? "n/a" : tokens[last_used]);
             }
 
             if (reactants) {
@@ -507,7 +506,11 @@ void setupReaction(Reaction& R, const AnyMap& node)
             reactants = false;
         }
     }
+}
 
+void setupReaction(Reaction& R, const AnyMap& node)
+{
+    parseReactionEquation(R, node["equation"]);
     // Non-stoichiometric reaction orders
     std::map<std::string, double> orders;
     if (node.hasKey("orders")) {
@@ -924,6 +927,44 @@ void setupElectrochemicalReaction(ElectrochemicalReaction& R,
     }
 }
 
+void setupElectrochemicalReaction(ElectrochemicalReaction& R,
+                                  const AnyMap& node, const Kinetics& kin)
+{
+    setupInterfaceReaction(R, node, kin);
+    R.beta = node.getDouble("beta", 0.5);
+    R.exchange_current_density_formulation = node.getBool(
+        "exchange-current-density-formulation", false);
+}
+
+bool isElectrochemicalReaction(Reaction& R, const Kinetics& kin)
+{
+    vector_fp e_counter(kin.nPhases(), 0.0);
+
+    // Find the number of electrons in the products for each phase
+    for (const auto& sp : R.products) {
+        size_t kkin = kin.kineticsSpeciesIndex(sp.first);
+        size_t i = kin.speciesPhaseIndex(kkin);
+        size_t kphase = kin.thermo(i).speciesIndex(sp.first);
+        e_counter[i] += sp.second * kin.thermo(i).charge(kphase);
+    }
+
+    // Subtract the number of electrons in the reactants for each phase
+    for (const auto& sp : R.reactants) {
+        size_t kkin = kin.kineticsSpeciesIndex(sp.first);
+        size_t i = kin.speciesPhaseIndex(kkin);
+        size_t kphase = kin.thermo(i).speciesIndex(sp.first);
+        e_counter[i] -= sp.second * kin.thermo(i).charge(kphase);
+    }
+
+    // If the electrons change phases then the reaction is electrochemical
+    for (double delta_e : e_counter) {
+        if (std::abs(delta_e) > 1e-4) {
+            return true;
+        }
+    }
+    return false;
+}
+
 shared_ptr<Reaction> newReaction(const XML_Node& rxn_node)
 {
     std::string type = toLowerCopy(rxn_node["type"]);
@@ -986,9 +1027,18 @@ unique_ptr<Reaction> newReaction(const AnyMap& node, const Kinetics& kin)
     }
 
     if (kin.thermo().nDim() < 3) {
-        unique_ptr<InterfaceReaction> R(new InterfaceReaction());
-        setupInterfaceReaction(*R, node, kin);
-        return unique_ptr<Reaction>(move(R));
+        // See if this is an electrochemical reaction
+        Reaction testReaction(0);
+        parseReactionEquation(testReaction, node["equation"]);
+        if (isElectrochemicalReaction(testReaction, kin)) {
+            unique_ptr<ElectrochemicalReaction> R(new ElectrochemicalReaction());
+            setupElectrochemicalReaction(*R, node, kin);
+            return unique_ptr<Reaction>(move(R));
+        } else {
+            unique_ptr<InterfaceReaction> R(new InterfaceReaction());
+            setupInterfaceReaction(*R, node, kin);
+            return unique_ptr<Reaction>(move(R));
+        }
     }
 
     if (type == "elementary") {

test/data/surface-phases.yaml

@@ -7,6 +7,7 @@ phases:
 - name: Pt-surf
   thermo: ideal-surface
   kinetics: surface
+  reactions: [Pt-reactions]
   species: [{Pt-surf-species: all}]
   site-density: 2.7063e-9 mol/cm^2
   state: {T: 900 K, P: 1 atm, coverages: {Pt(s): 0.5, H(s): 0.4, O(s): 0.1}}
@@ -17,6 +18,31 @@ phases:
   site-density: 5e-17 mol/cm
   state: {T: 1073.15 K, P: 1 atm, coverages: {(tpb): 1.0}}
 
+- name: graphite
+  thermo: lattice
+  species: [{graphite-anode-species: [C6, LiC6]}]
+  state: {T: 300, P: 101325, X: {C6: 1.0, LiC6: 1e-5}}
+  density: 2.26 g/cm^3
+
+- name: electrolyte
+  thermo: lattice
+  species: [{graphite-anode-species: [Li+(e), PF6-(e), EC(e), EMC(e)]}]
+  density: 1208.2 kg/m^3
+  state:
+    T: 300
+    P: 101325
+    X: {Li+(e): 0.08, PF6-(e): 0.08, EC(e): 0.28, EMC(e): 0.56}
+
+- name: anode-surface
+  thermo: ideal-surface
+  kinetics: surface
+  reactions: [graphite-anode-reactions]
+  species: [{graphite-anode-species: [(int)]}]
+  site-density: 1.0 mol/cm^2
+  state:
+    T: 300
+    P: 101325
+
 
 gas-species:
 - name: H2
@@ -82,7 +108,7 @@ tpb-species:
     model: constant-cp
 
 
-reactions:
+Pt-reactions:
 - equation: H2 + 2 Pt(s) => 2 H(s)
   rate-constant: {A: 4.4579e10 cm^3/mol/s, b: 0.5, Ea: 0.0}
   orders: {H2: 1.0, Pt(s): 1.0}
@@ -91,3 +117,43 @@ reactions:
   coverage-dependencies: {H(s): {a: 0, m: 0, E: -6000 J/mol}}
 - equation: H + Pt(s) => H(s)
   sticking-coefficient: [1.0, 0.0, 0.0]
+
+
+graphite-anode-species:
+- name: EC(e)
+  composition: {C: 3, H: 4, O: 3}
+  thermo: {model: constant-cp}
+- name: EMC(e)
+  composition: {C: 4, H: 8, O: 3}
+  thermo: {model: constant-cp}
+- name: O2(e)
+  composition: {O: 2}
+  thermo:
+    model: constant-cp
+    h0: 0.00177285 kJ/mol
+    s0: 143.5777069 J/mol/K
+- name: Li+(e)
+  composition: {Li: 1, E: -1}
+  thermo: {model: constant-cp}
+- name: PF6-(e)
+  composition: {P: 1, F: 6, E: 1}
+  thermo: {model: constant-cp}
+- name: C6
+  composition: {C: 6}
+  thermo: {model: constant-cp}
+- name: LiC6
+  composition: {C: 6, Li: 1, E: -1}
+  thermo:
+    model: constant-cp
+    h0: -11.65 kJ/mol
+# dummy species for anode/electrolyte interface
+- name: (int)
+  composition: {}
+  thermo: {model: constant-cp}
+
+
+graphite-anode-reactions:
+- units: {length: cm, quantity: mol}
+- equation: LiC6 <=> Li+(e) + C6
+  rate-constant: [5.74, 0.0, 0.0]
+  beta: 0.4

test/kinetics/kineticsFromYaml.cpp

@@ -221,3 +221,19 @@ TEST(Kinetics, InterfaceKineticsFromYaml)
     auto IR3 = std::dynamic_pointer_cast<InterfaceReaction>(R3);
     EXPECT_TRUE(IR3->is_sticking_coefficient);
 }
+
+TEST(Kinetics, ElectrochemFromYaml)
+{
+    shared_ptr<ThermoPhase> graphite(newPhase("surface-phases.yaml", "graphite"));
+    shared_ptr<ThermoPhase> electrolyte(newPhase("surface-phases.yaml", "electrolyte"));
+    shared_ptr<ThermoPhase> anode(newPhase("surface-phases.yaml", "anode-surface"));
+    std::vector<ThermoPhase*> phases{anode.get(), graphite.get(), electrolyte.get()};
+    auto kin = newKinetics(phases, "surface-phases.yaml", "anode-surface");
+    graphite->setElectricPotential(0.4);
+    vector_fp ropf(kin->nReactions()), ropr(kin->nReactions());
+    kin->getFwdRatesOfProgress(ropf.data());
+    kin->getRevRatesOfProgress(ropr.data());
+
+    EXPECT_NEAR(ropf[0], 0.279762523, 1e-8);
+    EXPECT_NEAR(ropr[0], 0.045559637, 1e-8);
+}