diff --git a/doc/sphinx/yaml/reactions.rst b/doc/sphinx/yaml/reactions.rst index 4af1a1d9a..fb7513160 100644 --- a/doc/sphinx/yaml/reactions.rst +++ b/doc/sphinx/yaml/reactions.rst @@ -28,8 +28,10 @@ The fields of a ``reaction`` entry are: - :ref:`Chebyshev ` Reactions on surfaces or edges are automatically treated as - :ref:`interface ` reactions, without the need to - specify the ``type``. + :ref:`interface ` reactions, and reactions which + involve charge transfer between phases are automatically treated as + :ref:`electrochemical ` 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 `_. + +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 diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index b0eaa9537..809089cba 100644 --- a/src/kinetics/Reaction.cpp +++ b/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 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 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 newReaction(const XML_Node& rxn_node) { std::string type = toLowerCopy(rxn_node["type"]); @@ -986,9 +1027,18 @@ unique_ptr newReaction(const AnyMap& node, const Kinetics& kin) } if (kin.thermo().nDim() < 3) { - unique_ptr R(new InterfaceReaction()); - setupInterfaceReaction(*R, node, kin); - return unique_ptr(move(R)); + // See if this is an electrochemical reaction + Reaction testReaction(0); + parseReactionEquation(testReaction, node["equation"]); + if (isElectrochemicalReaction(testReaction, kin)) { + unique_ptr R(new ElectrochemicalReaction()); + setupElectrochemicalReaction(*R, node, kin); + return unique_ptr(move(R)); + } else { + unique_ptr R(new InterfaceReaction()); + setupInterfaceReaction(*R, node, kin); + return unique_ptr(move(R)); + } } if (type == "elementary") { diff --git a/test/data/surface-phases.yaml b/test/data/surface-phases.yaml index 7295f627b..6f06e5070 100644 --- a/test/data/surface-phases.yaml +++ b/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 diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index 8d8259e7e..dfff7bf1c 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -221,3 +221,19 @@ TEST(Kinetics, InterfaceKineticsFromYaml) auto IR3 = std::dynamic_pointer_cast(R3); EXPECT_TRUE(IR3->is_sticking_coefficient); } + +TEST(Kinetics, ElectrochemFromYaml) +{ + shared_ptr graphite(newPhase("surface-phases.yaml", "graphite")); + shared_ptr electrolyte(newPhase("surface-phases.yaml", "electrolyte")); + shared_ptr anode(newPhase("surface-phases.yaml", "anode-surface")); + std::vector 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); +}