IdealGasPhase getStandardChemPotentials

[decaluwe]

System information

• Cantera version: 2.4
• OS: all (C++ source code)

Expected behavior The standard state definition for an ideal gas, typically, is pressure-independent. All pressure dependence in thermodynamic definitions is generally handled in the activity term, where the activity of a species k is equal to pressure*molefraction[k]/standardPressure.

More directly, this should always result in reverse reaction rate coefficients, calculated via k_rev = k_fwdexp(-deltaG^o/RT)prod(standardConcentrations^nu_k` being pressure independent (so long as f_fwd is pressure-independent).

Actual behavior (EDITED FROM ORIGINAL POST)

• For the IdealGasPhase, the standard state is defined as pressure-dependent (cf IdealGasPhase::getStandardChemPotentials, which includes muStar + log(pressure() / refPressure()) * RT();)
• For reactions where the forward rate k_fwd is supplied and the thermo is used to calculate k_ref (i.e. standard kinetics), this causes no problems. When it comes to calculating the reverse rate, the GasKinetics::updateKc routine includes the term m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc). In short, the logStandardConc term compensates for any pressure dependence in delta G^0.
• However, for interface reactions where rate_coeff_type = exchangecurrentdensity is used to supply the composition- and pressure-independent part of the exchange current density (let's call it i_o^*), which Cantera converts into a k_fwd and k_rev, we end up with reverse rate coefficients which are pressure dependent. @wbessler first noticed this when comparing charge-transfer reactions using forward and reverse rate constants (k_fwd and k_rev), versus using rate_coeff_type = exchangecurrentdensity. The former gives the expected behavior, whereas the latter does not.

Additional context See discussions below.

[decaluwe]

A few additional thoughts:

• There is a degree of guess as to exactly why the reverse rate is pressure-dependent for the electrochemical interface reactions, but I have not fully traced it through, yet.
• Regardless, I think the current approach in IdealGasPhase is flawed. I think overall the process ought to be more generalized, such that the chemical potential uses the phase's definition of activity for self-consistency, rather than hard-coding in the expression.
• I am happy to make these changes, but wanted to check in, first, to make sure my thinking on this is correct.

[bryanwweber]

I don't have enough knowledge to say if what you've got here is thermodynamically correct. We should get @cfgoldsmith and @rwest to weigh in, and maybe the SAB as well if you want.

That said, it would be interesting to see if you made these changes whether there would be any changes in the IdealGasPhase calculations. We expect not, but you never know... So I'd be interested in seeing an implementation 😁

[decaluwe]

There wouldn’t be any changes AFAIK— the bulk properties are all correct, it’s just how those are split between the “standard” and “non-standard” portions that causes problems for a heterogeneous reaction. Curious to see if it is a problem Generally, or if there is something specific about electrochemistry. For the Kc calculations, the changes to delta G^0 and logStandardConc would offset each other.

Okay, this has clarified that I’m pretty sure I’m right 😁, and I’ll at least work up an implementation to test it out.

Sent from my iPhone

On Apr 4, 2020, at 7:58 AM, Bryan W. Weber notifications@github.com wrote:

 I don't have enough knowledge to say if what you've got here is thermodynamically correct. We should get @cfgoldsmith and @rwest to weigh in, and maybe the SAB as well if you want.

That said, it would be interesting to see if you made these changes whether there would be any changes in the IdealGasPhase calculations. We expect not, but you never know... So I'd be interested in seeing an implementation 😁

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[speth]

I think it might be worth discussing the more general issue here, as we seem to keep bumping into this confusion about the differences between the standard / reference and unqualified thermodynamic properties. See, for instance https://github.com/Cantera/cantera/issues/522, where @decaluwe posed the question:

Do we need both (or either) "standard state" and "reference species property" calculations featured in the codebase? You can't really avoid having the reference data loaded and stored as such, but I'm somewhat skeptical that we need the standard state data explicitly calculated or stored.

I would go a step further, and question whether either "standard state" or "reference state" thermodynamic properties should be part of the public interface to the ThermoPhase class. These methods are used in a handful of places, such as the kinetics calculations mentioned above, as well as in the equilibrium solvers. However, it might be the case that these could be replaced with use of the full thermodynamic properties without losing anything.

[decaluwe]

1. @speth -- While I agree that this is a conversation/decision we keep putting off, and one that needs to occur, it is only partially relevant to the issue raised here (if I am understanding your suggestion correctly. I may not be - see below).

   Whether or not they are publicly accessible (I can think of 1 or 2 use cases where they are not necessary, but maybe a user would like them -- see below for this, too), the standard-state Gibbs energy of reaction is required internally for calculating the reverse reaction rate coefficient for reversible reactions. The issue here is that something seems to be awry with how this is done for certain reactions, such that Cantera calculates reverse rate coefficients that are pressure-dependent but ought not to be.

2. I appear to have misinterpreted @wbessler's original observation. The exchangecurrentdensity gives the incorrect (pressure-dependent k_rev) result, whereas the standard approach does not. I have amended the original issue message to reflect this. This makes sense (a bug of the nature I described above would have been caught by now, I would think), and makes the problem much less widespread (it does not impact all interface reactions, just this one rate coefficient type). It also suggests two possible paths to fixing it.

• Path 1: One could argue that the standard state for an IdealGasPhase should be pressure independent.
• Path 2: Another path is to simply look at how the exchangecurrentdensity is calculating k_rev. It likely assumes a pressure-independent standard state, and relaxing this assumption will likely fix things.

This actually breaks the issue into two separate issues:

1. Fixing the current issue: The conversion and calculation of k_rev is likely not coded in a generalized manner when rate_coeff_type = exchangecurrentdensity is used. I will inspect. Clearly, the code should just implement the full theory, without making any assumptions about the standard state, and let the terms be whatever they are.

   1a. On a related note, I did notice that something similar is done in IdealGasPhase::getChemPotentials - the routine hard-codes the species activities as X_k, rather than calling sub-routines to calculate them. It's done correctly, but if, for example, we do change the standard state, this is less robust in that it would not automatically adjust.

2. A broader question about the standard state for the IdealGasPhase, and whether it should be pressure-dependent. Pressure-independent would be more consistent with how the field conceives of the Ideal Gas phase, likely helping prevent similar confusion and or errors on two fronts:

• Users doing their own calculations involving the standard state properties and/or species activities might expect a constant P_o.

• Future developers might also assume that a pressure-independent standard state for this class (particularly since there is a separate class, for now, for IdealGasVPSS). Again, if they code in a generalized manner it shouldn't matter, but realistically there may be times where these errors sneak through (present case included).

  @speth, here is where your comment comes in: if we say "standard state properties are no longer publicly accessible", it fixes error type 1 (which is likely the more significant issue), with the slight inconvenience to any users who might want to access those properties (I would assume this is a small but non-zero number of users). It does not address the second type of error/confusion, and punts on the overall issue of what a developer might expect the phase to assume.

  Also, I might be mis-interpreting your suggestion,

  "However, it might be the case that these could be replaced with use of the full thermodynamic properties without losing anything."

  Are you suggesting that we would use the full Gibbs energy of reaction to calculate the reverse rate coefficient? FWIW, the standard state Gibbs energy of reaction is explicitly required in this calculation. If we eliminate the routine standard-state routines, it would require the k_rev routine back out the standard state Gibbs energy of reaction from the full term, as part of the calculation. Not impossible, but not terribly efficient and another possible source of error/confusion.

  These tradeoffs might be acceptable/preferable, but we ought to at least be fully aware of what the tradeoffs are.

[decaluwe]

So... I think this discussion on standard states is valuable, but is unrelated to the error in the charge transfer rate coefficients.

I went and re-derived the conversion, yesterday. The error observed by @wbessler turns out to be either a typo or simply an error in the derivation of InterfaceKinetics::convertExchangeCurrentDensityFormulation. There is a missing m_beta[i] exponent in the standard concentration term in the denominator.

I have a fix implemented and will submit a PR shortly. What to do with this discussion, though? I think it is worth figuring out what to do with the standard/reference/etc. states, even though it might not be at the very top of anyone's list.

[speth]

I think the enhancements repo is the place to start this discussion, and the fix for the issue you noted with convertExchangeCurrentDensityFormulation should close this issue. I can take care of creating the new issue in the enhancements repo. We can copy any relevant comments here or on other issues to that thread.