ck2yaml crashes when parsing a third body reaction with PLOG

[corykinney]

Problem description

ck2yaml crashes when parsing a third body reaction with PLOG.
I've tried searching open/closed issues for this, but either I'm searching with the wrong keywords or this is a new issue.

Steps to reproduce

1. Download AramcoMech 3.0 CHEMKIN files (zip file in supplementary materials)
2. Run ck2yaml --input="AramcoMech 3.0.MECH" --thermo="AramcoMech 3.0.THERM" --transport="AramcoMech 3.0.TRAN"
3. See error:

Error reading reaction starting on line 5544:
"""
C4H6+M<=>C2H3+C2H3+M                                         +1.3489600E+072 -1.6550000E+001 +1.4340000E+005
PLOG /                                       +3.9473700E-002 +2.3442300E+074 -1.7690000E+001 +1.4080000E+005 /
PLOG /                                       +7.8947400E-002 +2.2908700E+074 -1.7510000E+001 +1.4210000E+005 /
PLOG /                                       +1.5789500E-001 +4.3651600E+073 -1.7130000E+001 +1.4290000E+005 /
PLOG /                                       +3.1578900E-001 +1.3489600E+072 -1.6550000E+001 +1.4340000E+005 /
"""

ERROR: Unable to parse 'AramcoMech 3.0.MECH' near line 5544:
Reaction C4H6+M<=>C2H3+C2H3+M
 contains parameters for more than one reaction type.

System information

• Cantera version: 2.6.0a3
• OS: Windows 10
• Python version: 3.8.10

[ischoegl]

Hi @corykinney … thanks for reporting. To date, third-body collision partners are only supported for ThreeBodyReaction and FalloffReaction, which is why ck2yaml complains about this input.

[corykinney]

I was afraid that it might not be supported yet. Is there a recommended way to convert mechanisms with these reactions in them, or any plans to support them in the future, @ischoegl?

[ischoegl]

Not sure about a recommended way at the moment. From my perspective, adding this capability to 2.6 would make sense (e.g. most likely by allowing third-body collision partners for any reaction type supported by GasKinetics). This will likely require some discussion over at Cantera/enhancements as this would be a departure from past practice.

[bryanwweber]

@corykinney How does CHEMKIN handle calculating these rates? Does the rate calculation match what you'd expect from the equation, that is, kf*[C4H6]*P/RT ? I'm asking because I wonder why this hasn't been supported. It seems like something that would be fairly trivial to implement in the initial version of the PLOG calculation (in Cantera), so it seems odd to me that it would not be implemented unless this kind of reaction didn't make sense kinetically, and Chemkin's default behavior is to just ignore the +M part.

[speth]

When I implemented PLOG reactions in Cantera, I believe I checked what Chemkin did if a third body was specified and came to the conclusion that the rate constant is not multiplied by [M] regardless of how you write the reaction equation and that the effect of third bodies has to be folded into the rate constants. This is therefore what ck2cti does -- the reaction in question is converted to just:

pdep_arrhenius('C4H6 <=> C2H3 + C2H3',
               [(0.0394737, 'atm'), 2.344230e+74, -17.69, 140800.0],
               [(0.0789474, 'atm'), 2.290870e+74, -17.51, 142100.0],
               [(0.157895, 'atm'), 4.365160e+73, -17.13, 142900.0],
               [(0.315789, 'atm'), 1.348960e+72, -16.55, 143400.0])

But maybe it's worth checking again to see if Chemkin thinks C4H6 <=> C2H3 + C2H3 and C4H6 + M<=> C2H3 + C2H3 + M should be handled differently.

[ischoegl]

@speth ... interesting. If it's a quirk of CHEMKIN to ignore it, then ck2yaml should probably generate a warning but run through. The mechanism comes out of NUIG, so it's something Cantera should support without having to tweak. Funny enough it won't convert for me even with ck2cti, although a different error is thrown (and the mechanism dates to 2018). ... Edit: it does convert to the expression above, albeit only if transport is omitted.

[speth]

Sorry, I omitted the fact that I had to ignore transport to get the mechanism to convert with ck2cti -- I think that allowing things like 1.000 as the "geometry" flag is a little dodgy, but there are a handful of mechanism authors that generate input files looking like this.

[ischoegl]

Ok, great. AramcoMech 3.0 uses both PLOG with and without third-body collision partners, which by itself doesn’t necessarily imply anything. I don’t have a way to check what CHEMKIN does. If it ignores M as in the past, the ‘fix’ would be simple (I’d still advocate for a warning); if M is now considered by CHEMKIN, then ck2cti is no longer correct and would introduce artifacts. Implementing third body concentrations for other reaction types would not be difficult …

[corykinney]

@ischoegl If there are PLOG reactions both with and without third-body collision partners, to me that would imply, regardless of whether or not CHEMKIN ignores it, that the mechanism authors intend for it to be considered in the reaction. If it's not difficult to add third-body to other reaction types, I think that would be beneficial for improving the capabilities of Cantera, as long as it makes sense kinetically. If for some reason CHEMKIN does ignore it, perhaps there could be a flag on the ck2cti/ck2yaml scripts so the user can decide whether to replicate CHEMKIN's behavior by ignoring it or to model it as it is written.

[ischoegl]

@corykinney ... while it's not difficult to implement, the change also isn't trivial, so I believe that @speth's suggestion to first check what CHEMKIN does here is the best approach. Regarding using inconsistent conventions, AramcoMech 3.0 also uses the PLOG reaction HOCO<=>CO+OH, which would presumably involve a third-body collision partner also. My assumption here is that mechanism development uses an onion-peel approach, where new parts are added while core elements aren't necessarily updated.

[rwest]

There has been some discussion of issues like this in RMG circles in the past, eg.
ReactionMechanismGenerator/RMG-Py#1084 (comment)
ReactionMechanismGenerator/RMG-Py#1084 (comment)
so I'm tagging @alongd in case he recalls the outcome.
In my experience, guessing what mechanism authors intended is often hazardous. ¢2

[ischoegl]

@rwest ... thanks for the links! Regarding CHEB, specifying a third body collision partner is currently ignored by Cantera (with a warning being issued in the upcoming 2.6 release - see #1183). While I haven't worked on this, I am under the impression that this behavior is consistent with CHEMKIN?

For PLOG, the question on hand is what to do with the third-body collider in AramcoMech 3.0. @alongd: I'd be curious about the current CHEMKIN behavior with an unspecific collider (+M)?

Fwiw, there is a feature request that would allow for alternative third-body reaction types in Cantera/enhancements#133 (at the moment, they're only considered for ThreeBody and Falloff).

[alongd]

@ischoegl, I'm no longer a Chemkin user, but I can interpret my comment from a few years ago: Chemkin seemed to work fine when specifying a general collider (+M), but only allowed a specific collider, e.g., (+Ar), in CHEBYSHEV polynomials, but not in PLOG representations. I hope this helps, happy to discuss more

[speth]

Hi @alongd, thanks for weighing in. I'm wondering if you recall what the specific behavior was in any of these cases. That is, If a PLOG reaction is written as A+B(+M)=C+D(+M), is the forward rate of progress calculated by multiplying the effective third-body concentration or not? Do you have any idea what it was doing for A+B+M=C+D+M?

And I guess, similar questions for CHEB reactions.

[alongd]

Hi @speth, my understanding is that a PLOG reaction in Chemkin is written as A+B=C+D without a +M or a (+M), e.g., see here and here, and that the forward rate per P is not multiplied by [M], but instead it is already a function of P: ln(k) for a PLOG representation is interpolated as in Eq. 3-35 in the Chemkin manual.

CHEB reactions are expressed as A+B(+M)=C+D(+M), here's an example from an RMG output:

NO(5)+OH(9)(+M)=HNO2(75)(+M)                        1.000e+00 0.000     0.000    
    TCHEB/ 500.000   1500.000 /
    PCHEB/ 0.790     1.184    /
    CHEB/ 6 4/
    CHEB/ 2.879e+00    5.688e-02    -7.667e-04   -1.590e-15  /
    CHEB/ 3.364e+00    1.604e-02    2.211e-04    -4.482e-16  /
    CHEB/ -8.282e-02   -2.600e-04   5.993e-05    7.267e-18   /
    CHEB/ -2.974e-02   6.843e-04    3.490e-06    -1.912e-17  /
    CHEB/ -1.096e-02   3.678e-04    2.288e-06    -1.028e-17  /
    CHEB/ 4.528e-03    -4.952e-04   -1.532e-07   1.384e-17   /

This is a bit confusing since I also think that the interpolated rate is not multiplied by [M] either, since it is already expressed in terms of P as well as in Eq. 3-38 in the above manual (so I don't know why it is written differently in Chemkin or whether it "could" be written that way, and that's just how RMG does it; I do know it compiles well in Chemkin).

As for the difference between A+B(+M)=C+D(+M) and A+B+M=C+D+M, as @rwest once pointed to me I believe, the former means that at the high-pressure-limit the rate isn't multiplied by [M], and at the low pressure limit it does as in a Troe or a Lindemann expression. The later means that the rate coefficient is always multiplied by [M], as in a "ThirdBody" expression (some more on that here). The behaviour here is different than PLOG or CHEB since these expressions aren't explicit functions of P.

I hope I understood the questions correctly, happy to share anything else that I can.

[ischoegl]

While I am not a CHEMKIN user, I tried to investigate this issue with ANSYS 2021 R1 Chemkin (I only have access to the graphical interface). Starting with the example closed_homogeneous__ignition__delay, I swapped the mechanism and (eventually) got things to run using C4H6 as fuel. I ran two different versions:

• Original version of the mechanism linked on top
• Version with +M removed for the block of PLOG reactions involving C4H6

Comparing CSV files, I don't see any difference in output, indicating that the +M is an optional 'decoration'.

[ischoegl]

@speth ... based on this outcome (I repeated this test twice, just to make sure), I believe there now is the option to just strip the +M in ck2yaml, while providing a warning?

[speth]

Thanks for doing the dirty work of checking Chemkin's behavior in this case, @ischoegl. I think that sounds like a perfectly reasonable resolution, and I'm glad to know that we don't need to support any weird new reaction type.