Improve handling of electrochemical reactions

[speth]

Abstract

There are several ways of specifying the rate constant for electrochemical reactions which are used in practice. The aim of this enhancement is to support input of rate constants using those formulations that are convenient for users, while keeping the underlying implementation reasonably clean & efficient.

Motivation

Cantera currently has had an incomplete and not fully-implemented set of options for defining electrochemical reaction rates. The non-functional options (namely, those related to the Butler-Volmer formulation) were deprecated in Cantera 2.5 (See Cantera/cantera#749 and Cantera/cantera#750).

In the discussion of that issue, @decaluwe provided a summary of the theoretical bases for electrochemical reaction types, which can be found here: Charge.Transfer.in.Cantera.pdf.

Some additional context provided by @wbessler in the above-linked issue is also worth highlighting: https://github.com/Cantera/cantera/issues/749#issuecomment-559727736

Possible Solutions

To summarize the suggestions from both @decaluwe and @wbessler, it makes sense to support two general types of electrochemical reactions, one corresponding to Marcus theory / elementary reactions, and one for implementing different variants of the Butler-Volmer equation. In more detail:

Marcus theory:

• For the ‘standard’ treatment described in Sec. 2.1 of the above-linked PDF, things operate largely as they should.
• Users should be able to either pass in a k_fwd or pass in i0 and have it converted to an equivalent k_fwd
• Currently, an electrochemical reaction is identified because it either (a) contains a specification of beta, contains a specification of i0, or, in the absence of either of these, the reaction equation is used to identify a charge-transfer reaction.
• @ischoegl has suggested that a required flag should be added. @decaluwe says "I think this is a good idea, but should be optional, on the user end. Correctly evaluating the reaction driving forces is not a ‘special feature’ that a user needs to request, imho."
• @decaluwe also suggests: I’ve got half a mind that we can derive a method to fold electric potential effects directly into the species activities, such that electrochemical reactions don’t require any special handling. Here, electric potential effects would be part of the activity calculation, where this impact evaluates to multiplication by 1.0 for non-charge-transfer reactions. This requires a little more thought, as it has both theoretical and implementation hurdles. But it is worth exploring.

Butler-Volmer:

Inputs: The user should be able to:

• Pass in a k_fwd and have it converted to an equivalent i0.
• Pass in i0 (likely via Arrhenius parameters)
• Pass in i*0 for a known reference state (known species concentrations), along with a list of these reference concentrations [X*_k], which is then used to calculate i0.
• Pass in non-stoichiometric species orders for i0
• Set a flag that allows the user to specify use of the Tafel equation

Outputs: The user should also be able to query, for a given state:

• The exchange current density i0
• The overpotential eta.

@ischoegl has also raised concerns about the fact that checking whether a reaction involves charge transfer currently involves parsing the reaction equation from its string form twice. If desired, this duplicate parsing could be eliminated by making parsing of the reaction equation a separate step from instantiation of a Reaction object.

[decaluwe]

Just a brief update: we have started some work on this, and have working modules in place for three different charge transfer mechanisms: Butler-Volmer, Marcus theory, and Marcus-Hush-Chidsey.

They are in a pretty basic state, at the moment: concentration dependence is only partially implemented, and setting other properties like reaction orders is yet-to-be-implemented. But it is probably a good time to write up a basic test suite and submit an initial PR, before adding in the "next layer" of features.

[speth]

As a concrete example, some of this has been implemented in a branch (https://github.com/CSM-Offenburg/cantera/blob/a91ed333a0f6448f0ab99fa76c9be0f1368ce03a/src/kinetics/InterfaceKinetics.cpp#L349-L389). The following is a part of InterfaceKinetics::updateROP():

    for (size_t i = 0; i < m_beta.size(); i++) {
        // Get the reaction index, relative to the array of all interface 
        // reactions:
        size_t irxn = m_ctrxn[i];

        // Update electrochemical potentials. This is equal to the 
        //  overpotential, multiplied by the elementary charge transfered and 
        //  the Farday constant, for each reaction:
        vector_fp echemPotentials(nReactions(), 0.0);
        getDeltaElectrochemPotentials(&echemPotentials[0]);

        // Reciprocal of RT:
        double rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
        // Calculate the BV-type form directly, if specified:
        if (m_ctrxn_form[i] > 0 && m_ctrxn_form[i] < 4) {
            // Scale the concentration dependence by the refeence
            // concentrations:
            m_ropf[irxn] *= m_ctrxn_concScaleFactor[i];

            double beta = m_beta[i];
            if (m_ctrxn_form[i] == 3){ // marcus theory
                beta += echemPotentials[irxn]/4/m_lambda[i];
            } 
            // Calculate rop_f and rop_r:
            m_ropf[irxn] *= exp(-beta*echemPotentials[irxn]*rrt)/Faraday;
            m_ropr[irxn] = m_ropf[irxn]*exp(echemPotentials[irxn]*rrt);
        } else if (m_ctrxn_form[i] == 4) {
            // Scale the concentration dependence by the reference
            // concentrations:
            double i_o = m_ropf[irxn] * m_ctrxn_concScaleFactor[i];
            double eta_RT = echemPotentials[irxn]*rrt;
            double lambda_RT = m_lambda[i]*rrt;
            double erfc_term = erfc((lambda_RT - sqrt(1.0 + sqrt(lambda_RT)
                            + eta_RT*eta_RT))/2/sqrt(lambda_RT));
            m_ropf[irxn] = i_o*sqrt(Pi*lambda_RT)
                           * erfc_term/(1. + exp(eta_RT))/Faraday;
            m_ropr[irxn] = i_o*sqrt(Pi*lambda_RT)
                           * erfc_term/(1. + exp(-eta_RT))/Faraday;
        } 
    }

[decaluwe]

Okay, so moving discussion from #cantera1416 over here.

@korffdm has an implementation that works with the current ReactionRateFactory format. It is not the cleanest implementation, in terms of making the theory easy to trace and debug, but it is accurate (it passes local tests we have developed and will make part of the PR, if we stick with the current development), and appears to be the shortest development path to getting these capabilities into Cantera.

One note, w/r/t the theory described above.

• My earlier understanding of "Marcus theory" was somewhat confused. What I had previously described as Marcus theory was just an alternate implementation of Butler-Volmer theory. In short: Butler-Volmer assumes that the species energy scales linearly with the phase electric potentials. The two approaches described above both make this linear assumption, but implement it using a mass-action form and what we'll call a "Butler-volmer form." But they are both implementations of the same theory (as evidenced by the fact that they both reproduce the exact same rates for a given state).
• @korffdm's implementation implements the actual Marcus theory - where species energy is a quadratic function of the electric potential. This is in addition to adding Butler-Volmer (B-V) theory capabilities, where both Marcus and B-V use the B-V reaction format.

All that aside, I would like to discuss if an implementation closer to the code snippets in @speth's comment on Feb 15 is possible. From the refactorization, it looks as though all production rates currently go through updateROP, where the procedure is:

1. Update the rate coefficients for the current state (temperature, concentrations, electric potential effects)
2. Calculate reverse rate coefficient from k_fwd/K_eq
3. Multiply the fwd rate coefficient by the product of the activity concentrations raised to reactant species orders to obtain the forward ROP
4. Multiply the ref rate coefficient by the product of the activity concentrations raised to product species orders to obtain the reverse ROP

Before I speculate further, let me make sure I understand the situation correctly: am I correct that all reactions go through this procedure?