Mechanism validation fails for PLOG rates

[Cory Kinney]

I recently upgraded from Cantera 2.5.0 to Cantera 2.6.0, and the AramcoMech 3.0 mechanism no longer works with Cantera 2.6.0. The mechanism (available from NUI Galway) fails validation and generates the following error:

CanteraError thrown by PlogRate::validate:
Invalid rate coefficient for reaction 'C4H6 <=> C3H3 + CH3'
at P = 15999, T = 200.0
Invalid rate coefficient for reaction 'C4H6 <=> C3H3 + CH3'
at P = 31760, T = 200.0

The input for that reaction is:

C4H6+M<=>CH3+C3H3+M +4.0738000E+196 -5.0000000E+001 +2.4140000E+005
PLOG / +3.9473700E-002 +1.5849000E+148 -3.7240000E+001 +1.8850000E+005 /
PLOG / +7.8947400E-002 +8.9125000E+159 -4.0320000E+001 +2.0130000E+005 /
PLOG / +1.5789500E-001 +5.2481000E+196 -5.0000000E+001 +2.4360000E+005 /
PLOG / +3.1578900E-001 +4.0738000E+196 -5.0000000E+001 +2.4140000E+005 / 

And the resulting reaction after running ck2yaml is:

- equation: C4H6 <=> CH3 + C3H3 # Reaction 2176
type: pressure-dependent-Arrhenius
rate-constants:
- {P: 0.0394737 atm, A: 1.5849e+148, b: -37.24, Ea: 1.885e+05}
- {P: 0.0789474 atm, A: 8.9125e+159, b: -40.32, Ea: 2.013e+05}
- {P: 0.157895 atm, A: 5.2481e+100, b: -50.0, Ea: 2.436e+05}
- {P: 0.315789 atm, A: 4.0738e+100, b: -50.0, Ea: 2.414e+05}

I was going to make an issue on GitHub but figured I'd ask here first.

[Bryan Weber]

Hi Corey,

This may be a red herring, but I noticed that your A factors seem really high. If I recall correctly, 1e200 is almost certainly beyond the collision limit for a reaction rate. 

In any case, can you calculate the rate under those conditions by hand and see what you get?

Best,

Bryan
 

[Cory Kinney]

Bryan,

I can try to calculate it by hand to see what the rate looks like. I thought the same thing about the A factors; however, I tried changing the order of magnitude of those A factors to 1e100 and still received the same error, so it wouldn't appear to be the cause of the issue. I think it's also relevant that this is the AramcoMech 3.0 mechanism, which is fairly widely used in CHEMKIN with no issues. 

Cory

[Bryan Weber]

Hi Cory,

I think 1e100 is still much too high for the A-factor. My recollection (although I don't have a source) is that reaction rates above about 1e30 are above the collision limit, or at least are suspect. This paper lists an equation to calculate the collision limit: https://repository.kaust.edu.sa/bitstream/handle/10754/662466/Chatelain2020a_CLV_Preprint.pdf?sequence=1&isAllowed=y It does require some geometric information about the molecules in question, so I didn't dig any further into specific estimates. That paper, and I think reference [14] in that paper, also point out that there are dozens of reactions in many published chemical kinetic models whose reaction rates exceed the maximum rate that we believe is physically possible!

I also think that CHEMKIN's ability to load a model and run a simulation to the specified end time is only a weak indicator of the correctness of the result. The error that Cantera is reporting here (invalid rate constant) would also cause problems for CHEMKIN, except we choose to inform the user before their simulation fails. Or, if a CHEMKIN simulation does not fail, they must have implemented some undocumented check of their own to "fix" the rate constant under these conditions. For what it's worth, I suspect the reason that simulations run to completion with CHEMKIN is that the temperature is rarely set as low as 200 K; one can quibble with our choice to test reaction rates at that low of a temperature, but personally I don't see why reaction rates should be invalid at any point for a valid input value. This topic was also raised by Richard West, see: https://github.com/Cantera/enhancements/issues/54

Separately, I will again object strongly to using CHEMKIN as the benchmark for any calculation. I object on the principle of using proprietary and (more importantly) closed-source and publicly unverified code for open science research. And more specifically, I know of at least one bug in CHEMKIN related to handling negative A-factors in PLOG reactions, where CHEMKIN is able to either take the natural logarithm of a negative number, the equations in their own theory manual are not being implemented by the software, or some other normalization is used that they don't disclose.

I think the value of the validation checks that Cantera includes and which CHEMKIN apparently doesn't is self-evident from the fact that Cantera is telling you that the reaction rate calculated by these parameters is unphysical. I really think people using models and software need to be aware of the "garbage in/garbage out" principle... if you use a model with unphysical parameters, you should not expect physical results... and this is one way that Cantera does that. (Other ways are behind the --permissive flag in the input file converters, because we're also pragmatists in addition to idealists :-))

Best,

Bryan

[Cory Kinney]

Bryan,

I agree that CHEMKIN, as a closed-source proprietary software, shouldn't be the benchmark and Cantera shouldn't seek to replicate it. When I first found Cantera my mindset was that it was aiming to be an open-source version of CHEMKIN. Of course now I better understand the goals of Cantera; therefore, I agree with Cantera's philosophy to enforce physical constraints and raise errors when reaction rates are unrealistic to allow for better informed mechanism development. 

My dilemma is that, for better or for worse, reaction mechanisms are published with these unphysical rates, and even though they can't possibly be modeling the actual combustion chemistry accurately (as you said "garbage in garbage out"), that CHEMKIN output was fit to some dataset within reasonable uncertainty and I'd like to be able to at least get the predictions for a given set of conditions, even if the mechanism is inherently flawed. 

I'm sure its difficult balancing the goals of Cantera with maintaining some level of compatibility with CHEMKIN. Perhaps there could be an InvalidRateConstantWarning that could be raised when the rate constant at 200 K is invalid, so that the user is aware that the simulation might fail, and an InvalidRateConstantError that could be raised if the temperature were set to a value that would give an invalid rate constant? This way the simulation could still be run at a higher temperature where the rate might be valid but there would still be a warning that the reaction is physically incorrect?

Cory

[Bryan Weber]

Hi Cory,

Thanks for the reply! I'm glad you picked up that Cantera isn't an open-source CHEMKIN clone, I think that's really fundamental to the project and as you say, not totally obvious. If you have an opinion about how Cantera should handle this specific error, please feel free to comment on that issue I linked to and share your ideas :-)

Best,

Bryan